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CHAPTER  1 
INTRODUCTION 

With  the  advent  of  high-resolution  imaging  techniques,  it  has 
become  feasible  to  employ  pictures  taken  from  geostationary  satellites 
for  the  interpretation  and  prediction  of  atmospheric  processes  in  the 
mesoscale  range.  Although  polar-orbiting  satellites  with  a  sun- 
synchronous  orbit  have  a  relatively  low  flight  altitude  of  less  than 
1000  km  and  thus  provide  images  with  high  resolution,  areas  in  the  mid¬ 
latitudes  are  over-flown  only  twice  a  day.  Geostationary  meteorological 
satellites,  however,  like  the  American  GOES  or  the  European  Meteosat, 
orbit  in  36,000  km,  but  transmit  a  discretized,  digital  image  of  the 
full  earth  disc  each  half  an  hour,  with  subpoint  resolution  of  1  -  2.5 
km  in  the  visible  channel  and  5  -  8  km  in  the  thermal  infrared  (i.e., 
approximately  11  urn)  channel .  For  smaller  areas  of  coverage,  higher 
repetition  rates  can  be  selected  by  ground  control.  High  reception 
rates  of  up  to  every  3  minutes  allow  the  meteorological  image  analyst 
to  detect  fast  changing  features  that  bear  relevance  in  local  forecast¬ 
ing,  but  might  be  negligible  for  the  description  of  synoptic-scale 
situations. 

In  spite  of  the  introduction  of  satellite-based  radiometers  for 
other  spectral  bands,  such  as  the  water  vapor  channel  around  6  ym  on 
Meteosat,  or  microwave  sounding  devices  aboard  the  spacecrafts,  the 
classic  meteorological  analysis  from  satellites  is  based  on  the  inter¬ 
pretation  of  cloud  images.  Clouds  are  the  most  regularly  visible 
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phenomenon  providing  excellent  information  on  fundamental  conditions. 
The  monitoring  of  cloud  clusters,  vorticity  centers,  frontal  bands  as 
well  as  single  clouds  are  of  great  importance  to  the  local  forecaster. 
One  of  the  foremost  objectives  is  the  identification  and  prediction  of 
severe  storms.  For  a  detailed  forecast,  however,  it  is  essential  to 
obtain  the  cloud  parameters  not  only  qualitatively  but  also  quantita¬ 
tively. 

With  the  support  of  computerized  digital  image  processing 
facilities  the  parameters  of  interest  are  retrieved  and  displayed  such 
that  most  of  the  information  is  processed,  formatted  and  then  made 
available  to  the  analyst  who  operates  as  a  'man-in-the-loop. '  Most  of 
the  conclusions  are  based  on  the  motion  of  clouds  and  their  respective 
top  temperature.  Motion  and  size  parameters  of  a  cloud  can  be  employed 
to  estimate  wind  fields  [1,2]  as  well  as  to  establish  models  of  cloud 
development  and  extrapolate  from  past  observations. 

Current  cloud  displacement  estimators  comprise  the  subjective 
film  loop  technique  [3],  the  'correlation'  technique  [4,5]  and  the 
ISODATA  method  [6,7],  an  intricate  cluster  analysis  technique.  The 
objective  of  this  research  is  to  develop  estimators  that  reveal  not 
only  information  on  lateral  (x-y)  cloud  displacement,  but  also  on  ro¬ 
tation  and  size  changes.  The  latter  two  subjects  are  new  in  the  quan¬ 
titative  cloud  analysis.  Several  methods  will  be  introduced  in  this 
work,  based  on  different  signal  characterizations  and  various  degrees 
of  image  abstraction. 

In  Chapter  2  various  approaches  to  the  signal  representation 
of  cloud  images  are  discussed.  The  classical  method  of  estimating 
lateral  displacement  employs  the  luminance  (gray-tone)  information  of 
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the  discretized  digital  image.  The  estimate  is  obtained  by  identifying 
the  location  of  the  maximum  of  the  inner-product  surface  (for  energy 
signals)  or  of  the  spatial  cross-correlation  (for  power  signals)  of  two 
consecutive  pictures.  Improvement  of  the  maximum-peak  detectability  by 
means  of  prefiltering  is  discussed  in  Chapter  3.  In  Chapter  4  it  is 
shown  how  the  luminance  values  of  a  cloud  or  a  cloud  field  can  be 
fitted  to  a  bivariate  Gaussian  function.  The  lines  of  constant  density 
are  ellipses  about  the  center  of  luminance  gravity.  By  following  these 
ellipses  through  a  sequence  of  images,  clouds  may  be  traced  and  infor¬ 
mation  obtained  on  rotation,  displacement  and  size  changes.  Chapter  5 
deals  with  the  analysis  of  cloud  contours.  A  contour  is  traced  with 
constant  velocity.  In  equidistant  time  intervals,  the  values  x  and  y 
of  the  contour  points  are  combined  to  a  complex  number  z  =  x  +  jy. 

This  complex  series  can  be  transformed  into  a  Discrete  Fourier  Series 
(Fourier  descriptor).  The  resulting  coefficients  reveal  direct  infor¬ 
mation  on  lateral  and  rotational  displacement  as  well  as  the  scale 
change  of  the  contour.  A  numeric  matching  algorithm  determines  the 
amount  of  rotation,  shifting  and  scale  change  in  order  to  match  two 
contours  in  the  optimum  sense.  Matching  can  be  performed  even  on  the 
basis  of  a  truncated  Fourier  descriptor  series,  thus  reducing  the 
computational  requirements. 

The  inverse  transformation  of  a  descriptor  series  truncated  to 
the  three  lowest  coefficients  is  an  ellipse  approximating  the  original 
contour  in  the  least-squares  sense.  These  ellipses  can  be  tracked  and, 
if  desired,  be  compared  with  the  ellipses  obtained  from  the  bivariate 
Gaussian  luminance  fit. 
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The  problem  of  parameter  forecasting  is  approached  in  Chapter  6. 
A  Kalman  filter  is  introduced  to  track  cloud  parameters  and  obtain 
one-step  predicted  values,  such  as  the  geometric  parameters  of  an 
ellipse  from  a  given  series  of  elliptic  approximations.  Linear  and 
non-linear  models  and  the  respective  requirements  for  a-priori  infor¬ 
mation  are  discussed.  The  work  concludes  with  a  simple  case  study. 


CHAPTER  2 


SIGNAL  ANALYSIS 

2.1  The  Inner  Product 

A  digital  satellite  image  with  clouds  is  given  on  a  Cartesian 
matrix  grid  at  times  t  =  t . ,  i  =  0,  1,  2 ,  .  .  .by  the  image  sequence 
(x^n.m)},  n  =  0,  .  .  .  N  -  1;  m  =  0,  1  .  .  .  M  -  1,  (1) 

with  x  representing  the  non-negative  digitized  brightness  value,  e.g., 
an  integer  number  between  0  (=  black)  and  255  (=  white,  for  8-bit  rep¬ 
resentation).  For  displacement  analysis,  all  images  have  to  be  navi¬ 
gated  to  some  reference  mark,  and,  if  possible,  the  brightness  of  the 
clouds  adjusted  for  sun  angle  changes. 

The  image  (1)  is  a  combined  time-space  process  and  can  be  con¬ 
sidered  as  deterministic  in  time  and  space,  deterministic  in  space  but 
random  in  time  (or  vice-versa),  or  random  in  both  space  and  time.  Al¬ 
though  atmospheric  processes  seem  to  be  described  most  appropriately 
by  random  variables  and  stochastic  processes,  a  statistical  space 
characterization  of  the  cloud  image  (1)  is  problematic.  The  customary 
approach  of  defining  an  image  as  a  realization  of  a  two-dimensional 
Gaussian  random  field  is  often  too  simplistic  and  cannot  account  for 
the  existence  of  well  defined  edges,  contours  and  texture  [8].  Addi¬ 
tional  problems  arise  with  the  background.  It  cannot  be  considered  as 
additive  noise  since  it  is  partially  occluded  by  the  clouds.  This 
Implies  that  a  complete  statistical  description  of  the  time-space 
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process  (1)  cannot  be  easily  obtained.  Rather  than  regarding  each 
single  pixel  (i.e.,  picture  element)  of  (1)  as  belonging  to  the  domain 
of  the  process  realization,  the  following  view  has  been  adopted.  Each 
single  image  is  a  noisy  measurement  of  time-varying  parameters  such  as 
location,  size  and  cloud  top  temperature.  Instead  of  estimating  the 
statistics  of  each  pixel,  the  goal  is  to  (i)  extract  the  parameters  of 
interest  from  each  image  and  (ii)  develop  predictors  for  these  para¬ 
meters  on  the  basis  of  past  measurements. 

The  author  is  convinced  that  models  based  on  structural  para¬ 
meters  are  closer  to  the  aim  of  the  analysis  of  the  underlying  meteoro¬ 
logical  processes  and  therefore  preferred  to  the  classical  decomposition 
of  the  signal  into  orthogonal  functions. 

With  the  exception  of  snow-covered  terrain  and  occasional  sea- 
surface  reflection,  the  background  of  satellite  images  is  usually 
darker  than  the  clouds.  To  remove  disturbing  influences  on  displace¬ 
ment  analysis,  the  background  is  set  to  zero  (=  black)  by  either  sub¬ 
traction  or  clipping  at  a  threshold  level.  Snow-covered  ground  is 
identified  by  comparison  with  the  reflected  radiance  in  the  infrared 
(IR)  bands  [9].  Furthermore  let  us  assume  that  the  cloud  of  interest 
is  completely  defined  within  the  image  frame.  Since  the  cloud  image 
sequence  (x(n,m)}  then  is  zero  outside  some  area  of  support  0<n<N-l, 
0<m<M-l,  it  belongs  to  the  class  of  finite  area  sequences  [10].  Finite 
area  sequences  are  energy  signals,  i.e., 

0  <  |j  x(n,m)  |i  2  <  ®  (2) 

with  the  energy  of  x(n,m)  defined  as 

2  N=1  M=1  2 

|(x(n,m)||  =  <  x(n,m),  x{n,m)  >  =  L  l  x£(n,m) 

n=o  m=o 


(3) 
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The  operation  ||  .  .  .  ||  also  is  called  the  norm,  and  <.  .  .>  denotes 
the  inner  product.  These  signals  describe  the  ^  space.  The  inner 
product  of  two  sequences  {x(n,m)}  and  {y(n,m)}  with  both  areas  of 
support  given  as  O^nfN-l,  0<m^M-l,  is  defined  as 
Kxy{0,0)  =  <x{n,m),  y(n,m)> 

M-l  N-l 

=  l  l  x(n,m)  y(n,m).  (4) 

m=0  n=0 

If  the  sequence  {y(n,m)}  is  displaced  by  v  and  u,  we  obtain 
M-l  N-l 

K  (v,u)  *  1  1  x(n,m)  y(n-v,m-u) .  (5) 

xy  m=0  n=0 

If  v,p  in  (5)  are  varied  from  0.  .  .N-l,  0.  .  .M-l,  there  is  at  least 

one  set  of  indices  v  and  u  such  that  K  is  maximum.  For  the  case  that 

xy 

(y(n,m)}  is  a  shifted  replica  of  (x(n,m)},  i.e., 

(y(n,m) }  =  {x{n-mo,m-mo)}t 

the  absolute  maximum  of  <xy(.)  is  obtained  for  0  »  nQ,  v  =  mQ. 

Hence  the  inner  product  K  (.)  can  be  used  to  identify  the  same  cloud 

xy 

in  two  consecutive  images  and  estimate  its  lateral  displacement. 

Although  K  (.)  is  commonly  placed  among  cross-correlation  esti- 
xy 

mators,  the  inner  product  is  not  based  on  any  probabilistic  descriptions. 
The  cross-correlation  of  a  random  field  is  defined  as  a  statistical 
expected  value  type  of  operation  involving  the  joint  probability  density 
function  [11].  However  a  random  field  of  energy  signals  is  dynamic 
enough  that  statistical  Inferences  may  not  be  made  by  the  analysis  of  a 
single  image.  Displacement  estimation  with  the  inner  product  (5)  is  a 
classical  operation  in  image  processing  [12-14].  In  spite  of  the  absence 
of  any  statistical  assumptions  it  is  commonly  called  the  "cross¬ 
correlation"  method. 
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Sometimes  the  displacement  of  features  is  to  be  estimated  within 
scenes  that  cannot  be  described  as  energy  signals  but  as  a  space- 
truncated  power  signal.  Assume  that  the  feature  image  {x(n,m)>  of 
size  L,  K  is  to  be  matched  within  a  large  search  image  {y(n,m)}  of 
size  N,  M.  In  those  instances  a  maximum  of  the  inner  product  (5)  would 
not  necessarily  indicate  optimum  fit.  Hence  the  normalized  'correlation- 
coefficient1  surface  [4,  J*  12-15]  is  introduced  as 


L-l 

K-l 

Z 

Z  x(l  , 

k)  yO-v,  k-u) 

1=0 

k=0 

•L-l 

K-l  , 

lh  TL-l  K-l  ?  1*5 

Z 

Z  x*(l 

,k)  Z  Z  y^(l-v.k-y) I 

ll=0 

rr 

II 

o 

J  Li=o  k=0  J 

(6) 


The  performance  of  this  operator  depends  on  the  size  of  the  feature 
image  and  the  local  statistics  of  both  feature  and  search  area  [16]. 

If  the  sequences  (x(n,m)}  and  (y(n,m)}  are  finite  area  sequences  with 
identical  extensions,  the  denominator  of  (6)  is  a  constant,  and  the 
numerator  is  equal  to  the  inner  product  (5). 

Since  the  background  removal  is  often  desired,  the  influence  of 
this  operation  on  the  inner  product  needs  careful  consideration.  Back¬ 
ground  removal  by  clipping  is  performed  according  to 
<x(n,m)  for  x(n,m)  <  c 

x  (n,m)  =  (7) 

'0  otherwise 

with  c  a  specified  threshold.  Background  is  subtracted  by 
ix(n,m)  -  c  for  x(n,m)  <  c 


xs(n,m) 


0  otherwise 


(8) 


Since 

xc(n,m)  =  x$(n,m)  +  c, 

the  inner  products  of  {xc(n,m)}  and  (xs(n,m)}  are  related  as 
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<xc(n,m),  xc(n,m)>  =  <xs(n,m),  xs(n,m)> 


+  <c-q(xs(n,m)),  c-q(x$(n,m))> 


with  the  plateau  function 


q(xs(n,m))  = 


1  for  x$( . )  >  0 


(9) 


(10) 


'0  otherwise  . 

Equation  (9)  implies  that  the  inner  product  of  the  clipped  signal  is 
equal  to  the  inner  product  of  the  subtracted  signal  plus  the  inner  prod¬ 
uct  of  the  plateau  function  c-q(.)  and  hence  introduces  a  bias  term 
depending  on  the  threshold  level.  Therefore  the  subtraction  method 
with  subsequent  dynamic  range  widening  is  preferred. 

In  spite  of  the  author's  perseverance  in  the  distinction  of 
inner  product  and  cross-correlation,  a  particular  random  model  can  be 


assumed  for  finite  area  sequences  that  establishes  identity  of  both 
operations.  With  the  finite  area  sequence  {x(n,m)>  being  zero  outside 
the  rectangle  0<n<N-l,  0<m<M-l,  the  following  random  process  is  con¬ 
structed: 

X(n,m)  =  X(n-V ,m-W)  (11) 

with  V,  W  two  independent  random  variables  of  uniform  distribution, 
i.e.,  V  =  U[0,Zn]  and  W  =  U[0,Zm].  The  parameters  Zp  and  Zm  define 
the  limits  of  the  random  displacement  of  the  sequence  {x(n,m)}.  The 
statistical  cross-correlation  of  this  process  is  given  by 
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N  M 


VN  VM 


l  Z  x,  (n,m)  x.  (n-f;,  m-n))/  Z  Z  q(x,  (n,m) -x,  (n-£,m-n) ) , 
n=0  m=0  1  n=0  m=0  1 


with  the  q-function  as  defined  in  (10). 

The  expected  value  of  this  expression  then  is,  for  the  given  probability 
distribution,  the  cross-correlation 
£{X(n,m)  X(n-5,m-n)} 

N  M 

=  1  Z  Z  x(n,m)  x(n-£,m-n).  (14) 

(Zn+N)(zm+M)  n=0  m=0 

This  result  is  identical  to  the  inner  product  (5)  times  a  constant 
factor.  Remember  that  this  relationship  holds  only  for  the  random 
displacement  with  uniform  distribution.  Since  the  motion  of  clouds  as 
well  as  many  other  physical  objects  is  not  governed  by  a  uniform  dis¬ 
placement  probability  but  rather  by  dynamic  equations,  the  short-comings 
of  the  inner  product  as  correlation  estimator  are  apparent. 


2.2  Spectrum  and  Moments 


For  finite  area  sequences,  the  two-dimensional  discrete  Fourier 
transform  (DFT)  is  defined  as  [10,  p.  117] 

X(p,q)  =  DFT  {x(n,m)> 

=  Z  Z  x(n,m)  exp  {-2ttj(£^  +  ^J)}  (15) 

n=0  m=0  1 

where  j  =  /T,  and  the  variables  p,q  are  the  spatial  frequency  wave 
numbers.  The  numbers  X(p,q)  are  the  Fourier  series  coefficients  of  a 
periodic  extension  of  (x(n,m)}.  From  the  definition  (15)  it  follows 
that  the  two-dimensional  DFT  is  a  separable  transform,  i.e.,  it  can  be 
implemented  as  a  sequence  of  one-dimensional  DFT's.  The  sequence 
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{X ( p ,q ) >  is  also  referred  to  as  the  Fourier  spectrum,  or  briefly, 
spectrum  of  (x(n,m)>.  It  can  be  written  as 

X(p,q)  *  Ax(p,q)  exp  {-J>x(p,q)}  (16) 

with  Ax(.)  the  amplitude  and  4>x(p,q)  =  arg  [X(p,q)]  its  phase. 
Alternatively,  X(p,q)  can  be  expressed  as 

X(p«q)  =  XR(p,q)  +  j  Xj (p,q) ,  (17) 

XR(.)  denoting  the  real  part  and  Xj(.)  the  imaginary  part.  Amplitude 
and  phase  is  commonly  replaced  by  the  magnitude 

|X(p,q)  |  =  (XR2(p,q)  +  X^p.q))*4  (18) 

and  the  principle  phase 

_i  X. (p,q) 

ARG  [X(p,q)]  =  tan  ^  (19) 

R 

Note  however  that  magnitude  and  principle  phase  are  ambiguous  with  re¬ 
spect  to  amplitude  and  phase:  (i)  the  magnitude  is  non-negative  by 
definition  and  thus  introduces  phase  discontinuities  of  +it  whenever  the 
amplitude  changes  its  sign;  (ii)  the  inverse  tangent  function  leads  to 
phase  discontinuities  of  +ir  or  +2n  according  to  its  definition  *). 

The  product  of  the  spectra  of  two  finite  area  sequences  (x(n,m)} 

and  (y(n,m)}  is  the  Fourier  transform  of  the  inner  product  K  (v,u). 

xy 

If  {y(n,m)>  is  a  replica  of  (x(n,m)},  but  shifted  by  nQ  and  mQ,  the 
cross-spectrum  becomes 


*)  Most  FORTRAN  libraries  contain  two  different  inverse  tangent 
functions:  ATAN(X)  and  ATAN2(A,B).  Our  definition  of  tan-*  coincides 
with  the  two-argument  function  ATAN2(A,B).  The  difference  is  illus¬ 
trated  in  the  following  table: 


l  tt/4 

tt/2 

3tt/4 

TT 

5tt/4 

3tt/2 

7ir/4 

2"  i 
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X(p,q)Y*(p,q)  =  A^(p,q)  exp  +  -§•)}  (20) 

Hence  the  displacement  estimation  could  be  performed  equivalently  in 
the  frequency  domain  by  evaluating  the  phase  of  the  cross-spectrum. 

For  reasons  stated  above,  the  principle  phase  is  wrapped  around  multi¬ 
ple  values  of  n  and  therefore  does  not  permit  immediate  visualization 
of  the  phase  behavior.  The  evaluation  of  the  continuous  phase,  also 
called  phase  unwrapping,  has  gained  considerable  interest  in  the 
numerical  evaluation  of  the  cepstrum  [10,  p.  508;  17,  18].  The  cepstrum 
is  defined  as  the  inverse  Fourier  transform  of  the  complex  logarithm  of 
the  spectrum.  The  existence  of  two-dimensional  cepstra  is  shown  in 
[19].  In  the  appendix  of  this  work  a  two-dimensional  phase-unwrapping 
technique  is  introduced  based  on  Tribolet's  method  [17]  for  one¬ 
dimensional  signals.  The  phase  of  the  cross-spectrum  is  related  to  the 
spatial  displacement  of  spectral  energy.  If  both  input  images  differ 
only  by  a  lateral  displacement,  the  phase  is  linear  (see  (20)).  How¬ 
ever  if  the  images  display  other  changes  besides  displacement,  such  as 
size  and  shape  changes,  the  phase  terms  of  the  input  spectra  do  not 
differ  just  by  a  linear  component. 

There  are  some  interesting  relations  linking  inner  product, 
spectrum  and  moments  with  each  other.  Defining  the  j-k-  moment  of 
(x(n,m)}  as 

N-l  M-l  •  . 

s •  l  l  nJ  m  x(n,m) ,  (21) 

J  n=0  m=0 

then  the  total  brightness  or  luminance  mass  of  the  image  is  given  by 
Sqq.  The  center  of  luminance  gravity  is  derived  as 

.  S01  _  s10 

Pm  Sqq  '  Pn  Sqq 


(22) 
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The  partial  derivatives  of  a  two-dimensional  spectrum  are  evaluated  as 
^^1=  DFT  {-jnx(n.m)}  =  X'Rp(p.q)  +  jX'p(p,q) 

(23) 

dc)  =  DFT  =  ^  Rq(P*q)  +  Jxfq<P*q) . 

with  the  subscripts  R  for  real  part,  I  for  imaginary  part,  p  as 
derivative  with  respect  to  p,  and  q  with  respect  to  q.  The  derivatives 
of  phase  and  amplitude  then  can  be  derived  as 

d  arg  X(p,q)  _  XRXIp  *  XRpXI  .  . 

dp  x2  +  x2 

*R  *1 


dp 


X  X'  + 

Vru 


X  X" 

Viu 


(X2R  + 


x2)4 


(25) 


(similarly  with  respect  to  the  variable  q).  A  numerical  integration 
of  (24)  leads  to  an  approximation  of  the  unwrapped  phase.  With  the 
spectrum  at  zero  frequencies  being  purely  real,  i.e., 

N-l  M-l 

X(0,0)  ■  £  E  x(n,m)  =  XR(0,0)  =  sQQ  (26) 

n=0  m=0 

and  its  derivative 


d  X(0,0)  _ 
dp 


N-l  M-l 

-j  E  E  nx(n,m) 
n=0  m=0 


-jXfu(0,0) 


(27) 


the  derivative  of  the  phase  at  zero  frequency  is  given  by 

d  arg  X(0,0)  _  XlV0,0)  _  ^10 

dp  X^OTO]  s00  ’ 


(28) 


i.e.,  the  phase  derivative  with  respect  to  p  is  identical  to  the  nega¬ 
tive  value  of  the  location  of  the  luminance  center  of  gravity  in  n- 
direction.  The  same  holds  for  the  derivative  with  respect  to  q,  being 
equal  to  the  negative  value  of  the  location  in  m-direction.  Similar 
results  can  be  derived  for  the  phase  of  the  cross-spectrum.  Here  the 
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derivative  of  the  phase  at  zero  frequency  is  identical  to  the  negative 
values  of  the  displacement,  i.e.,  the  difference  of  the  position  of 
the  centers  of  luminance  gravity. 

As  we  had  seen  before,  the  phase  of  the  cross-spectrum  is  linear 
for  images  being  identical  except  for  lateral  displacement.  The  deriva¬ 
tives  then  are  constant  for  all  frequencies. 

It  is  obvious  that  the  displacement  of  the  center  of  luminance 
gravity  is  very  often  a  satisfactory  indication  of  moving  objects.  As 
we  have  shown  in  this  section,  this  displacement  can  be  evaluated  in 
the  space  or  frequency  domain.  The  assumption  of  scene  changes  con¬ 
sisting  solely  of  lateral  displacement,  however,  has  to  be  judged 
carefully  for  its  validity  in  the  case  of  application  being  under 
research. 


CHAPTER  3 


IMPROVEMENT  OF  THE  INNER-PRODUCT  ME1 HOD 
3.1  Problems  in  the  Estimation  of  Cloud  Motion 

The  inner  product  method  as  presented  in  the  previous  chapter 
associates  the  amount  of  lateral  displacement  with  the  location  of  the 
absolute  maximum  of  the  inner  product  of  two  consecutive  images.  Cloud 
size  and  shape  changes,  cloud  rotation,  non-uniformity  of  wind  vectors 
in  different  altitudes,  to  name  only  few  of  the  typical  phenomena  the 
image  analyst  has  to  cope  with,  establish  serious  limitations  to  the 
performance  of  this  estimator  [20  -  22].  Whereas  in  more  advantageous 
situations  the  absolute  value  of  the  inner  product  maximum  is  a  clear 
indicator  for  having  matched  the  same  cloud  in  two  consecutive  frames, 
the  analysis  of  more  complex  scenes  might  require  a  separate  matching 
decision  prior  to  the  displacement  estimation,  either  by  the  human  oper¬ 
ator  or  a  very  sophisticated  automated  computer  analysis. 

The  inner  product  might  not  produce  meaningful  results,  if  the 
cloud  has  undergone  severe  changes.  In  those  instances  we  refer  to 
methods  introduced  in  Chapter  4  and  5.  Special  consideration  should  be 
given  to  the  rotation  of  a  cloud,  since  it  often  indicates  severe 
weather.  Rotation  therefore  is  not  regarded  as  a  factor  of  disturbance 
but  as  a  parameter  to  be  estimated  quantitatively.  In  these  instances 
it  is  suggested  to  estimate  the  rotation  first,  then  realign  both 
images  and  employ  lateral  displacement  estimation  on  the  rotation- 
corrected  images. 
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Rotation  could  be  estimated  by  an  inner  product  based  on  a  polar- 
coordinate  grid  point  system  centered  at  the  axis  of  rotation.  Since 
the  resulting  inner  product  would  be  very  sensitive  to  the  choice  of 
this  axis  point,  a  combined  lateral  and  rotational  inner  product  has 
to  be  evaluated  until  the  absolute  maximum  is  found.  Notwithstanding 
the  problems  of  interpolation  between  the  two  systems  of  coordinates, 
such  strategy  requires  excessive  computer  power.  One  of  the  require¬ 
ments  given  for  this  work,  however,  was  that  the  developed  algorithms 
were  feasible  for  mini-computers  with  memory  sizes  of  64  to  128K  byte. 
This  restriction  encroaches  already  on  the  performance  of  the  Cartesian 
inner  product.  The  maximum  image  size  available  with  the  equipment 
used  for  this  work  was  only  128  by  128  pixels,  and  even  this  could  be 
achieved  only  by  a  complicated  I/O  overhead  structure  and  using  the 
image  memory  for  the  storage  of  intermediate  results.  In  the  light  of 
the  very  satisfactory  performance  of  the  rotation  estimators  in  Chapter 
4  and  5,  further  consideration  is  directed  now  solely  towards  the 
improvement  of  the  lateral  displacement  inner  product  method. 

3.2  Pre-Filtering  by  Inverse  Filters 

As  we  have  indicated  in  the  last  section,  the  analysis  of  cloud 
motion  has  to  cope  with  the  fact  that  the  objects  of  interest  undergo 
changes  be- ween  two  consecutive  image  frames.  For  the  inner  product 
estimator  therefore  the  question  arises,  how  the  location  of  the  maximum 
of  the  inner  product  surface  is  related  to  tne  characteristics  of  the 
input  signals.  In  other  words,  how  reliable  is  the  inner  product 
method  given  the  fact  that  mapping  the  first  into  the  second  image  is 
inhomogeneous?  Certainly,  the  more  points  are  mapped  by  a  homogeneous 
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vector  field,  i.e.,  with  parallel  direction  and  equal  absolute  value, 
the  more  energy  is  concentrated  around  the  maximum  peak  of  the  inner 
product.  As  a  conclusion,  two  performance  goals  of  the  inner  product 
method  are: 

(i)  the  displacement  of  only  meaningful  features  should  be 
estimated, 

(ii)  the  maximum  peak  of  the  inner  product  surface  should  be  as 
pronounced  as  possible,  i.e.,  concentrate  as  much  energy 
as  possible. 

The  following  paragraphs  demonstrate  how  both  tasks  can  be 
achieved  simultaneously  by  a  single  technique. 

A  well-known  property  of  the  discrete  Fourier  transform  is  that, 
if  a  signal  is  not  too  severely  contaminated  by  aliasing  (i.e.,  by 
undersampling),  a  'narrow*  inner  product  peak  is  equivalent  to  a 
'broad'  cross-spectrum,  and  vice-versa.  The  limiting  case  is  defined 
by  the  inverse  transform  of  a  constant  (i.e.,  'white')  spectrum 
X(p,q)  =  1,  for  all  p,q.  The  inverse  transform  of  this  spectrum 
results  in  an  inner  product  with  a  single  impulse  at  v  =  0,  y  =  0,  i.e., 

K(v,y)  =  6(u,u) 

with  6(.)  being  the  unit  impulse  function 
(  1  for  v  =  0,  u  =  0 

6(v,u)  =  J  (29) 

'0  otherwise. 

The  inverse  transform  of  a  spectrum 

n  m 

X(p,q)  =  exp  (- j27T(— +  -^)) 
is  the  shifted  6- function 

1  for  v  *  n  ,  u  =  m 

6{v  -  nQ,  u  -  mQ)  = 


0  otherwise. 


(30) 


From  a  comparison  with  (20)  and  (16)  it  appears  that  (30)  is  the 
inner  product  of  identical  signals,  each  with  a  white  spectrum,  shifted 


laterally  by  n^  and  m  . 

o  o 

Since  our  images  of  interest,  in  general,  do  not  possess  a  white 
spectrum,  we  conclude  that  a  narrow  inner  product  peak  can  be  obtained 
only  if  the  input  signals  are  passed  through  a  'whitening'  filter 
prior  to  the  inner  product  operation  (Fig.  1). 


(x(n,m)}  * 

(y(n,m)}  o- 

Input  Prefilter  Inner  Product 

Figure  1.  Inner  product  of  pre-filtered  images. 

A  two-dimensional  discrete  filter  is  a  mapping  L  which  maps  an 
input  image  {x(n,m)>  into  an  output  image  (y(n,m)}  =  {L[x(n,m)]>.  The 
response  of  the  filter  to  the  6-function  (29)  is  called  the  impulse 
response  or  spread  function  of  the  filter 
h(n,m)  =  L[6(n,m)]. 

The  output  of  a  linear  shift-invariant  (LSI)  filter  is  a  discrete  con¬ 
volution  of  input  and  spread  function,  i.e., 

OO  00 

y(n,m)  -  r.  I  x ( k ,  1 )  h(n-k.m-l) 

k=-oo  l=-oo 
oo  00 

=  i  i  x(n-k,m-l)  h(k,l).  (31) 

k=-co  1  —  — oo 

The  two-dimensional  z-transform  of  an  array  of  complex  numbers  f(n,m) 
is  defined  as  [23] 


tu 
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F(z1,z2)  =  X  z  f(n,m)  (32) 

n=-«>  m=-°° 

where  z^  and  z2  are  complex  variables.  The  z-transform  evaluated  at 
the  surface  z^  =  expCjw^),  z2  =  exp(jw2),  is  equal  to  the  (continuous) 
Fourier  transform  of  {f (n ,ni) }  provided  this  array  is  absolutely  summable. 
If  (f(n,m)}  is  a  finite  area  sequence,  the  discrete  Fourier  transform 
(15)  is  identical  to  the  Fourier  transform  evaluated  at  the  frequen¬ 
cies  =  2rrp/N  and  w2  =  2irq/M. 

The  z-transform  of  the  output  of  a  digital  LSI  filter  is  given 
by  the  product  of  the  z-transforms  of  input  and  spread  function,  i.e., 
Y(z1,z2)  =  X(z1,z2)  H(z1,z2).  (33) 

The  z-transform  of  the  spread  function,  i.e.,  H(zj,z2)  is  called 
the  system  function.  From  (33)  follows  that  the  desired  system  function 
of  the  inverse  filter  should  be 


H(zrz2)  -  N-1  M-1 

I  Z  x(n,m)  z,  z9 
n=0  m=0 


(34) 


provided  the  inverse  exists. 

A  pre-filter  designed  to  improve  the  detectibil ity  of  the  inner 
product  maximum  peak  therefore  depends  on  the  spectrum  of  the  input 
signal.  For  one-dimensional  signals,  this  approach  has  been  suggested 
by  using  various  numerical  inversion  schemes  based  on  the  knowledge  of 
the  input  spectrum  [24].  One  disadvantage  of  these  techniques  is  that 
the  spectrum  of  the  signal  has  to  be  known  a-priori.  Furthermore,  the 
pre-filter  has  to  be  re-calculated  whenever  the  spectrum  of  the  input 
changes. 

In  the  following  section  it  will  be  shown  how  a  two-dimensional 
discrete  filter  can  be  designed  with  an  extended  Wiener  technique, 


approximating  the  inverse  filter  to  any  desired  degree.  We  shall  see 
that  a  low-order  Wiener  filter  will  be  a  satisfactory  approximation  for 
a  common  class  of  images,  and  hence  its  numerical  values  will  not  have 
to  be  changed  that  often.  Another  feature  of  this  filter  is  computa¬ 
tional  efficiency.  However  before  introducing  the  two-dimensional 
Wiener  filter  design  technique,  some  points  about  the  existence  of  such 
an  inverse  filter  are  in  order. 

Assume  the  inverse  filter  is  causal,  i.e.,  its  output  depends 
only  on  the  'present'  and  'past'  values.  This  terminology  is  of  course 
a  convention  used  in  the  theory  of  one-dimensional  time  functions.  In 
this  context,  'past,'  'present'  and  ’future’  are  defined  as  orienta¬ 
tions  in  space:  past  values  of  (x(n,m)}  are  all  values  (x(k,l)}  with 
k<n,  l<m;  future  is  defined  by  k>n,  l>m.  Call  { hc (k ,1 ) } ,  k>0,  1>_0, 
the  spread  function  of  this  inverse,  ca.isa1  filter.  Since  the  system 
function  (34),  in  general,  is  a  two-dimensional  oolynomial  of  infinite 
order,  the  spread  function  will  be  of  infinite  duration: 

H(z,,z?)  ~  Z  Z  h  (k,l)  z,‘k  z/1.  (35) 

1  £  k=0  1=0  c  1  c 

The  filter  (35)  is  called  stable  [25]  if  and  only  if 

CO  CO 

Z  Z  |h  (k,l ) |<  «.  (36) 

k=0  1=0  L 

Since  the  two-dimensional  z-transform  of  any  array  f(m,n)  converges  if 
and  only  if 

no  co 

Z  Z  |f(m,n)  z^  m  Zg  n|  <  °°»  (37) 

m=-oo  n=-«= 

we  conclude  that  the  LSI  filter  is  stable  if  and  only  if  its  system 
function  converges  for  jz^j  =  jzgj  =  1. 
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In  one -dimensional  systems  analysis,  the  proof  of  stability 
generally  is  established  if  the  poles  lie  inside  the  unit  circle.  From 
the  inversion  (34)  it  is  obvious  that  this  is  the  case  whenever  the 
zeros  of  the  signal  x(.)  lie  inside  the  unit  circle,  i.e.,  x(.)  is 
minimum-phase.  For  one-dimensional  functions,  the  location  of  the 
zeros  with  respect  to  the  unit  circle  decides  whether  the  signal  is 
minimum-phase,  maximum-phase  or  mixed-phase.  Hence  for  a  one-dimensional 
minimum-phase  signal  there  is  a  stable,  causal  inverse  filter  [26].  If 
the  signal  is  maximum-phase,  i.e.,  all  zeros  are  outside  the  unit- 
circle,  a  'left-going'  z-transform  and  an  anti-causal  system  function 
can  be  defined  to  assure  stability  and  convergence. 

However  a  two-dimensional  polynomial,  in  general,  cannot  be 
factored  into  a  product  of  first  order  polynomials  of  a  single  variable, 
i.e.,  there  is  no  equivalent  decomposition  in  poles  and  zeros.  The 
fundamental  theorem  of  algebra  does  not  exist  for  two-dimensional  func¬ 
tions  (with  the  exception  of  separable  functions).  The  principle  of 
minimum-phase  therefore  has  to  be  modified  for  two  dimensions  [27]. 

Definition.  A  sequence  (x(k,l)}  will  be  called  minimum-phase 
if  the  following  conditions  are  satisfied: 

(1)  The  two-dimensional  z-transform  Xfz^.Zg)  =  0,  evaluated  at 
any  jzjj  =  1,  has  no  zeros  outside  the  Z2  unit  circle. 

(2)  The  two-dimensional  z-transform  X(zj,z,,)  =  0,  evaluated  at 
any  | |  =  1.  has  no  zeros  outside  the  z\  unit  circle. 

Stability  and  minimum-phase  then  are  connected  by  the  following 


two  theorems: 
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Theorem  1  (Shanks;  [27  -  28]) 

Given  that  X(zj,z2)  is  a  polynomial  in  (z^Zg),  for  the  coef¬ 
ficients  of  expansion  of  1/X(z^,z2)  in  negative  powers  of  z^  and  z 2 
to  converge  absolutely,  it  is  necessary  and  sufficient  that  X(zj,z2) 
not  be  zero  for  |z  J  and  |z2|  simultaneously  greater  than  one. 

Theorem  2  (Shanks  [27]) 

A  sequence  {x(n,m)}  is  two-dimensional ly  minimum-phase  if  and 
only  if  the  polynomial  X(z^,z2)  =  0  has  no  zeros  for  which  | z -^ | >_1  and 
j  Z2 1^.1  simultaneously. 

From  Theorems  1  and  2  it  can  be  concluded  that  the  inverse  fil¬ 
ter  H(z^,z2)  =  1/X(z^,z2)  is  stable  if  X(zj,z2)  is  minimum-phase.  If 
X(z1>z2)  is  not  minimum-phase,  then  H(Zj,z2)  will  not  be  stable. 

However  if  X(z^,z2)  is  maximum-phase,  there  is  a  stable  inverse 

filter 

H(z,  ,z? )  =  E  l  h  ( k ,  1 )  z.K  z  '  (38) 

1  2  k=Q  1=0  a  12 

with  h  (k,l)  the  anti-causal  or  anticipatory  impulse  response.  The 

a 

output  of  this  filter  depends  on  present  and  future  input  values. 
Whereas  in  the  analysis  of  time-varying  functions  future  values  are 
not  always  available,  anti-causal  system  functions  do  not  lead  to  any 
problem  in  image  processing  as  long  as  the  complete  image  is  available 
at  the  time  of  processing. 

Since  cloud  images  are  generally  mixed  phase  type  of  signals, 
the  spread  function  of  a  stable  inverse  filter  will  be  both  causal  and 
anti-causal ,  i .e. , 

00  -k  -i 

H(z, »z?)  =  I  l  h ( k ,1 )  z,  z?  1 
1  c  k=-oo  !  =  -»  1  i 


(39) 


The  output  sequence  therefore  will  be  given  by  a  double-sided 
convolution 


y(n,m)  =  £  £  h(k,l)  x(n-k,m-l).  (40) 

k=-oo  1  =-co 

Note  that  the  input  sequence  is  of  finite  area,  the  double  sum  there¬ 
fore  will  actually  have  finite  boundaries. 

In  general  it  will  not  be  possible  to  find  a  finite  (i.e., 
realizable)  filter  that  exactly  inverts  the  spectrum  of  the  input  sig¬ 
nal.  Our  goal  therefore  consists  of  finding  a  filter  that  approximates 
the  inversion  according  to  some  criterion  of  optimality.  Our  choice  of 
filter-type  is  guided  by  the  convolution  (40)  and  its  advantages  of 
computational  implementation.  Furthermore  we  are  guaranteed  to  obtain 
a  stable  filter.  The  input-output  relationship  of  such  a  finite  impulse 
response  (FIR)  filter  is 
K  L 

y(n,m)  =  I  £  a(k,l)  x(n-k,m-l)  (41) 

k=-K  1=-L 

with  a(.)  the  causal/anti -causal  spread  function.  Note  that  we  expect 
the  values  of  a(k,l)  be  different  from  the  low-order  terms  of  the 
spread  function  (h(k,l)}  of  the  exact  inversion.  In  the  next  section 
it  is  shown  how  the  desired  filter  coefficients  (a(k,l)},  -K<k<K, 

-L<1<1,  can  be  derived  by  a  least-squares  minimization  technique. 

3.3  Design  of  the  Planar  Least  Squares  Inverse  Filter 

In  general,  it  will  not  be  possible  to  find  a  filter  (41)  that 
produces  exactly  the  desired  unit  pulse  image  with  a  white  spectrum. 

The  objective  of  this  section  is  to  find  filter  coefficients  (a(k,l)} 
-K<k<K,  -L<1<L,  such  that  its  output  approximates  the  unit  pulse  array 
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<5(n,m)  in  the  least  squares  sense.  This  approximation  method  was 
developed  by  Gauss  and  successfully  applied  in  one-dimensional  signal 
processing  [29].  In  the  sequel  it  will  be  shown  how  this  technique  can 
be  extended  to  a  two-dimensional  ( i . e . ,  planar)  least  squares  inverse 
(PLSI)  filter.  Since  the  least-squares  design  technique  is  often  con¬ 
fused  with  the  so-called  Wiener  filter  design  method,  it  should  be 
noted  that  the  former  is  not  based  on  any  statistical  assumptions, 
whereas  the  latter  requires  the  estimation  or  a-priori  knowledge  of 
spatial  correlation  functions  [30,  p.  228]. 

The  PLSI  filter  can  be  thought  of  as  a  cancellation  or  predic¬ 
tion  error  filter  [31,  32] 

K  L 

y(n,m)  =  x(n,m)  -  Z  Z  b(k,l)  x(n-k,m-l)  (42) 

k=-K  1=-L 
(k,l)r»(0,0) 

with  the  prediction  filter  coefficients  b(k,l).  Comparing  this  nota¬ 
tion  with  the  filter  (41),  one  can  see  that 
a(0,0)  =  1 

a ( k ,  1 )  =  -b(k,l )  for  (k,l)  t  (0,0).  (43) 

Recall  that  we  require  y(n,m)  =  6(n,m).  Hence  it  is  reasonable 
to  require  y(0,0)  =  x(0,0)  and  y(k,l)  for  (k,l)  f  (0,0)  as  small  as 
possible.  With  the  total  squared  error  thus  being 

0°  o° 

E  =  Z  Z  (y(n,m) )  ,  (44) 

n=l  m=l 

the  filter  coefficients  are  obtained  by  minimizing  E  with  respect  to 
each  coefficient  a ( k ,1 ) .  With  (42)  and  (43)  the  error  (44)  becomes 

«  »  /  K  L  \2 

E  =  Z  Z  (  Z  Z  a(k,l)  x(n-k,m-l)j 

n*l  m=l  \k=-K  1=-L  / 
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K  L  K  L 

=  E  E  E  E  a*(k,l)  a(h,1)  • 

k=-K  1=-L  h=-K  i=-L 

oo  oo 

E  E  x(n-k,m-l)  x(n-h,m-i),  (45) 

n=l  m=l 

with  *  denoting  the  complex  conjugate.  Since  the  double  sum 

OO  OO 

E  E  x{n-k,m-l)  x(n-h,m-i) 
n=l  m=l 

is  the  inner  product  Kxx(h-k,i-l ) ,  (45)  can  be  written  as 

K  L 

E  =  a*(-K,-L)  E  E  a ( h , i )  K  (h+K,i+L)  +  .  .  . 
h=-K  i=-L  xx 

K  L 

+  a*(0,0)  E  E  a ( h , i )  K  (h,i)  +  .  .  . 

h=-K  i=-L  xx 

K  L 

+  a*(K,L)  E  E  a(h.i)  K  (h-K.i-L).  (46) 

h=-K  i=-L  xx 

Setting  the  derivatives  of  (46)  with  respect  to  the  filter  coefficients 

a(f,g)  to  zero,  we  obtain  the  system  of  equations 

K  L  -K<f<K 

E  E  a ( k ,  1 )  K  { k-f ,  1  -g )  =  0  for  -Ug<L  .(47) 

k=-K  i=-l  xx  but  Tf.g)  t  (0,0) 

Inserting  (47)  into  (46),  the  minimal  error  becomes 

E»i"  ■  *  ■  JL  iL a<Mi  K|ui(k,,)'  <48) 

Equations  (47)  and  (48)  are  the  two-dimensional  normal  equations.  In 
matrix  form  they  can  be  written  as 
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with  R.  being  the  Toeplitz  sub-matrix 

J 

U  (DJI- -  KXX(-KJ) 

•  • 

■  • 

•  ■ 

•  • 

•  • 

•  • 

Rj  -  »WM) - Kxx(0,j) 

•  « 

«  • 

■  ■ 

•  ■ 

•  a 

. lSc,<^»' 

the  sub- vectors  of  filter  coefficients 


and  the  vector  q  is  a  column  vector  containing  zero  elements  except  for 
the  minimal  error  e  in  the  position  corresponding  to  a{0,0): 


~Kxx(-2K,j) 

a 

• 

a 

a 

• 

a 

Kxx(-K,j)  (50) 

a 

a 

a 

a 

a 
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Since  the  filter  coefficients  have  to  be  normalized  such  that  a(0,0)  =  1, 
the  value  e  in  (52)  can  be  replaced  by  a  1.  The  solution  for  the  optimum 
filter  coefficients  is  simply  achieved  by  taking  the  middle  column  of 
the  inverse  of  the  R  -  matrix  in  (49)  and  equating  it  to  the  filter 
coefficients  of  the  a-vector  in  (49).  Since  the  R-matrix  of  (49)  is 
block-Toeplitz,  its  inversion  can  be  done  recursively  by  a  modified 
Levinson  algorithm  [57].  The  predictor  filter  coefficients  { b( k , 1 ) > 
of  (42)  then  are  obtained  as 

b(j  ,k)  =  -a(j,k)/a(0,0).  (53) 

As  a  special  case  consider  the  class  of  images  with  an  exponential  inner 
product  of  the  form 

Kxx(n,m)  =  Kxx(0,0)  (54) 

with  p  a  constant  such  that  0  <  p  <  1. 

These  images  are  usually  referred  to  as  first-order  Markov 
processes,  although  a  strict  definition  would  require  the  replacement 
of  the  inner  product  by  the  spatial  correlation  function  [30,  p.  108]. 

The  constant  p  is  the  'correlation'  of  adjacent  pixels  in  rows  and 
columns. 

Real-world  images  approximate  the  Markov  property  quite  fre¬ 
quently  [33].  The  inner  product  sequence  is  a  Toeplitz  matrix  and 
bears  the  advantage  of  easy  diagonal ization  and  approximation  to  a 
circulant  matrix.  A  whitening  filter  for  such  a  Markovian  image  then 
is  a  function  of  the  single  parameter  p.  With  the  technique  of  solving 
the  normal  equations  (49),  a  filter  of  order  3x3  results  in  the 
matrix  of  filter  coefficients 
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1 


a (1,1)  a(l,0)  a(l,-l)  ' 

P2  -p(l+p2)  p2 

Ca(j,k)]  = 

a(0,l)  a(0,0)  a(0,-l) 

= 

-p(1+p2)  (1+P2)2  -p(1+p2) 

a(-l,l)  a(-l,0)  a(-l,-l) 

p2  -p(l+p2)  p2 

J 

(55) 


For  p  =  1,  i.e.,  high  adjacent  'correlation,'  the  filter  (55)  reduces 
to  the  Laplacian  edge  enhancement  filter  [11,  p.  482] 


[a(j,k)]  = 


f  1  -2  1 

2  4-2 

1  -2  1 


(56) 


l  i 

At  this  point  some  comments  about  the  connection  between  inverse 
filtering  and  edge  detection  are  in  order.  Since  many  actual  images 
display  a  high  amount  of  energy  in  the  lower  spatial  frequencies,  as 
long  as  the  signal -to-noise  ratio  is  sufficiently  high,  we  expect  an 
inverse  filter  to  have  high-pass  characteristics.  However,  the  high- 
frequency  energy  of  images  is  usually  due  to  sudden  changes  of  the 
brightness  value,  as  in  edges  etc. 

The  inner  product  between  two  different  cloud  images  that  are 
processed  through  such  a  filter  then  is  based  more  on  the  high-frequency 
energy  of  the  input  signals,  proportional  to  the  order  of  the  filter. 

The  notion  of  lateral  displacement  estimation  based  on  the  high- 
frequency  information  had  been  discovered  earlier  and  leads  to  the 
suggestion  of  passing  the  images  through  high-pass  filters  prior  to  the 
inner  product  operation  [14,  20,  34].  The  technique  described  in  this 
section,  however,  adapts  itself  to  the  images  under  consideration  and 
leads  to  a  filter  in  the  least  squares  optimum  sense. 

The  importance  of  edges,  contours  and  texture  as  reference  marks 
for  the  estimation  of  cloud  displacement  has  been  st3,  !  before  [4]. 

i 
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As  we  shall  see  in  the  next  section,  the  inner  product  of  a  pre-filtered 
set  of  images  leads  to  a  better  maximum  peak  performance,  as  long  as 
the  displacement  of  the  high-frequency  image  contents  is  rather  homo¬ 
geneous.  If  this  is  not  the  case,  i.e.,  if  no  part  of  the  high-frequency 
information  can  be  matched  with  any  other  part  of  the  second  input  frame, 
the  method  described  above  might  not  lead  to  any  significant  maximum  at 
all.  In  many  cases  of  practical  interest,  however,  the  change  of  cloud 
shape  occurs  only  in  some  parts.  The  inner  product  of  pre-filtered 
images  then  would  indicate  the  displacement  of  those  high-frequency 
parts  that  have  translated  in  a  more  homogeneous  fashion. 

3.4  Results 

To  demonstrate  the  effect  of  inverse  pre-filtering  on  the  inner 
product  surface,  three  pairs  of  consecutive  cloud  images  were  chosen 
with  various  degrees  of  changes  (Fig.  2(a-f)).  The  cloud  pair  I  of 
Fig.  2(a-b)  shows  almost  no  changes  besides  lateral  translation.  Cloud 
pair  II  (Fig.  2(c-d)  reveals  substantial  size  and  shape  changes  in  the 
North-East  section  of  the  cloud.  Image  pair  III  of  Fig.  2(e-f)  con¬ 
sists  of  a  field  of  loosely  connected  clouds  with  a  rather  fuzzy 
appearance. 

The  inner  product  surfaces  were  obtained  as  follows.  Upon  in¬ 
spection  of  the  histogram  of  the  given  images,  a  threshold  was  selected 
and  the  background  removed  by  subtraction.  It  should  be  noted  that  in 
the  case  of  infrared  (IR)  imagery,  the  removal  of  dark  background  areas 
is  equivalent  to  the  specification  of  a  temperature  (and  thus,  in  gen¬ 
eral,  an  altitude)  level.  Subsequent  estimation  of  lateral  cloud  dis¬ 
placement  then  would  be  based  on  the  information  corresponding  to 
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Fig.  2  (a,b).  Cloud  pair  I. 


Fig.  2  (c,d).  Cloud  pair  II. 


Fig.  2  (e,f).  Cloud  pair  III. 
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temperature  lower  than  the  given  threshold,  i.e.,  for  higher 
clouds. 

A  3  x  3  inverse  filter  then  was  derived  by  solving  the  normal 
equations  (49).  As  indicated  in  the  preceding  section,  the  filter 
coefficients  were  normalized  such  that  a(0,0)  =  1. 

The  input  images  were  pre-filtered  by  a  direct  implementation 
of  the  convolution  (41).  Subsequently  the  inner  product  was  calculated 
by  the  fast  Fourier  transform  (FFT)  method.  Since  the  equipment 
available  for  this  operation  consisted  of  a  PDP  11/60  mini -computer  and 
a  single-frame  Ramtek  9050  display  processor,  the  size  of  in-core 
memory  was  severely  limited.  Therefore  a  two-dimensional  FFT  program 
was  written  that  makes  use  of  the  Ramtek  image  memory  for  the  storage 
of  intermediate  floating-point  values  and  thus  avoids  time-consuming 
I/O  operations  with  mountable  storage  devices.  With  this  configuration 
a  maximum  image  size  for  the  inner  product  of  128  x  128  was  achieved. 

Since  the  process  of  pre-filtering  might  result  in  negative 
pixel  values,  a  constant  bias  was  added  and  later,  during  the  inner 
product  operation,  subtracted  again. 

In  Fig.  3(a-c)  the  inner  product  surfaces  are  presented  without 
prior  inverse  filtering;  Fig.  3  a  corresponds  to  the  image  pair  I,  Fig. 
3  b  to  pair  II  and  Fig.  3  c  to  image  pair  III.  The  smoothness  of  the 
surface  is  clearly  discernible.  In  spite  of  the  fact  that  the  location 
of  the  maximum  peak  can  be  evaluated  numerically,  its  interpretation  as 
a  meaningful  displacement  estimate,  however,  remains  questionable. 

The  inner  product  surfaces  in  Fig,  4  (a-c)  are  the  equivalent 
results  when  pre-filtering  with  a  3  x  3  PLSI  filter  is  performed  prior 
to  the  inner  product  operation.  In  Table  1  the  respective  (auto)  inner 
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Table  1.  Inner  product  surfaces  Kxx  (v,p). 


Image:  Fig.  2  (a) 


0.2287 

0.3712 

0.5451 

0.7037 

0.7883 

0.4096 

0.6014 

0.7875 

0.8755 

0.8006 

0.6325 

0.8319 

1.0000 

0.8319 

0.6325 

0.8006 

0.8755 

0.7875 

0.6014 

0.4096 

0.7883 

0.7037 

0.5451 

0.3712 

0.2287 

Image:  Fig.  2  (c) 

0.5319 

0.5943 

0.6504 

0.6720 

0.6582 

0.6269 

0.7335 

0.8201 

0.8085 

0.7336 

0.7237 

0.8771 

1.0000 

0.8771 

0.7237 

0.7336 

0.8085 

0.8201 

0.7335 

0.6269 

0.6582 

0.6720 

0.6504 

0.5943 

0.5319 

Image:  Fiq.  2  (e) 

0.5309 

0.6001 

0.6559 

0.6502 

0.6086 

0.6617 

0.7755 

0.8490 

0.8012 

0.7021 

0.7470 

0.9025 

1.0000 

0.9025 

0.7470 

0.7021 

0.8012 

0.8490 

0.7755 

0.6617 

0.6086 

0.6502 

0.6559 

0.6001 

0.5309 

product  surface  Kxx(v,p)  around  the  center  (v,p)  =  (0,0)  is  shown.  The 
values  are  normalized  with  respect  to  Kxx(0,Q)  such  that  the  affinity 
of  the  given  images  to  the  Markov  processes  (54)  can  be  established. 

The  3  x  3  inverse  filter  coefficients  are  listed  in  Table  2. 

Although  the  maximum  peak  of  Fig.  4  (a)  is  much  more  pronounced 
than  in  Fig.  3(a)  for  cloud  pair  I,  the  same  cannot  be  said  for  cloud 
pair  II.  As  Fig.  4(b)  shows,  the  method  of  pre-filtering  can  result  in 
a  multitude  of  relative  maxima.  Although  the  location  of  the  absolute 
maximum  is  still  in  the  vicinity  of  the  maximum  of  Fig.  3(b),  such  a 
result  does  not  necessarily  increase  our  confidence  in  the  practicability 
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Table  2.  PLSI  filter  coefficients. 


a ( 1 , 1 )  a ( 1 ,0 ) 

a ( 1  »-l )  ' 

a(0,l)  a(0,0) 

a(G,-l) 

a ( - 1 , 1 )  a(-l,0) 

a(-l ,— 1 ) 

4 

Imaqe:  Fiq. 

2(a) 

0.223 

-.471 

0.086 

-.338 

-.338 

0.086 

-.471 

0.223 

Image:  Fiq. 

2(c). 

0.136 

-.247 

0.046 

-.444 

-.444 

0.046 

-.247 

0.136 

Imaqe:  Fig. 

2(el 

0.249 

-.440 

0.137 

-.445 

1.000 

-.445 

0.137 

-.440 

0.248 

as  an  estimator  for  lateral  displacement.  The  link  between  high-pass 
filtering,  inverse  filtering  and  edge  enhancement  was  established  in 
the  previous  section.  The  result  of  Fig.  4(b)  then  can  be  interpreted 
as  follows.  Since  the  cloud  has  undergone  a  significant  change  in  some 
areas,  the  inner  product  of  the  high-frequency  segments  of  the  signal 
(such  as  edges  and  texture)  indicate  relative  match  at  several  vectors 
of  displacement.  In  other  words,  if  the  displacement  of  the  high- 
frequency  parts  is  not  homogeneous,  the  number  of  relative  maxima  will 
increase.  A  similar  result  can  be  seen  in  Fig.  4(c)  for  cloud  pair 
III.  Although  the  displacement  has  occurred  in  a  more  homogeneous 
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fashion  (resulting  in  a  more  pronounced  absolute  maximum),  the  quali¬ 
tative  effect  is  the  same. 

As  a  solution  to  the  foregoing  problem,  the  following  modifica¬ 
tion  of  the  filter  coefficients  is  suggested.  The  coefficients  a ( k , 1 ) 
for  all  indices  besides  (0,0)  are  pre-multipl ied  by  a  factor  c  with 
0<c<l.  This  factor  determines  the  amount  of  the  original  image  added 
to  the  pre-filtered  output.  In  essence  the  frequency  characteristic 
of  the  filter  changes  from  high-pass  to  high-emphasis.  For  c  =  0.5, 
the  filter  output  will  consist  of  the  normalized  sum  of  the  filtered 
image  and  the  original  input.  In  other  words,  the  choice  of  the  empha¬ 
sis  factor  c  determines  how  much  energy  of  the  lower  frequency  range 
is  to  be  subjected  to  the  inner  product  operation. 

In  Fig.  5(a-c)  an  emphasis  factor  of  0.75  was  chosen.  This 
allows  the  additional  leakage  of  25%  of  the  lower  frequencies  into  the 
maximum  peak  of  the  inner  product  surface.  The  results  therefore  rep¬ 
resent  a  compromise  between  maximum  peak  sharpness  and  detectibil ity. 

The  numerical  location  of  the  absolute  maximum  of  the  inner 
product  surface  plots  is  listed  in  Table  3.  Since  each  maximum  falls 
on  a  discrete  grid  point,  a  quadratic  surface  is  fit  to  the  discrete 
inner  product  values  around  the  maximum  and  a  search  is  done  to  find 
the  location  of  the  maximum  of  the  continuous  quadratic  surface. 

Since  the  grid  points  of  the  satellite  images  under  study  are 
spaced  with  an  equivalent  surface  distance  of  5  .  .  .  20  km,  a  high 
accuracy  in  the  determination  of  the  maximum  is  desired.  For  example 
it  is  desired  to  estimate  wind  speeds  from  cloud  displacement  with  an 
accuracy  of  +  3  m/sec. 
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Table  3.  Displacement  estimates  including  quadratic  interpolation 
(in  pixel  units). 


Cloud 

Pair 

I 

II 

III 

■  “  - 

An 

Am 

An 

Am 

An 

Am 

no  pre-filter 

-2.0 

1.2 

2.5 

0.6 

2.8 

4.1 

PLSI  pre-filter 

-1.9 

1.2 

-3.1 

2.0 

4.1 

4.0 

PLSI  pre-filter 
emphasis  factor 
c  =  0.75 

with 

-1.8 

1.1 

1.75 

0.8 

4.2 

3.7 

3.5  Summary 

This  chapter  dealt  with  the  improvement  of  the  maximum  peak 
detectibi 1 ity  of  the  inner  product  as  lateral  displacement  estimator 
for  clouds.  Since  a  narrow  inner  product  peak  is  equivalent  to  a 
broad  cross-spectrum,  pre-filtering  with  an  inverse  filter  is  suggested 
prior  to  the  inner  product  operation.  For  a  filter  order  much  lower 
than  the  image  size,  only  an  approximation  of  the  inverse  can  be 
achieved.  The  design  of  the  (two-dimensional)  planar  least  squares 
inverse  (PLSI)  filter  is  introduced.  Since  clouds  are  usually  mixed- 
phase  signals,  the  optimum  filter  is  not  causal.  This  does  not  lead 
to  any  problems  in  digital  image  processing.  The  filter  coefficients 
are  obtained  by  solving  the  normal  equations  extended  for  two  dimen¬ 
sions.  In  general,  the  pre-filter  obtains  high-pass  characteristic  and 
thus  emphasizes  edges  and  the  texture  of  cloud.  It  can  be  shown  that 
a  3  x  3  filter  designed  for  an  image  with  exponential  inner  product  is 
equal  to  an  edge  detector  found  in  the  i  iterature. 

The  method  of  PLSI  pre-filtering  is  demonstrated  by  3  examples 


of  consecutive  cloud  images.  As  long  as  the  pattern  of  the  high- 
frequency  information  does  not  change  drastically,  very  good  results 
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can  be  achieved.  Otherwise  it  is  suggested  to  modify  the  filter 
coefficients  by  multiplication  with  a  constant  factor.  This  changes 
the  frequency  characteristic  of  the  filter  from  high-pass  to  high- 
emphasis.  The  choice  of  the  factor  determines  the  amount  of  lower 
frequencies  exercising  influence  on  the  build-up  of  the  maximum  peak. 
Thus  a  compromise  between  peak  detectibil ity  and  sharpness  can  be 
achieved.  To  further  improve  the  identification  of  the  location  of  the 
maximum  peak,  the  inner  product  around  the  obtained  discrete  maximum 
value  is  fit  to  a  quadratic  surface,  and  its  maximum  searched  for. 


CHAPTER  4 

BIVARIATE  NORMAL  FIT 

4.1  Introduction 

From  the  results  of  the  foregoing  chapter  it  becomes  apparent 
that  the  inner  product  method  for  lateral  displacement  estimation  does 
not  perform  satisfactorily,  if  the  cloud  scene  has  undergone  severe 
changes  of  its  shape.  The  same  problem  arises  when  the  cloud  rotates, 
either  by  horizontal  shear  or  due  to  a  curved  track.  Since  the  inner 
product  in  the  form  presented  in  this  work  is  based  on  a  summation  of 
cross-energy  along  the  Cartesian  system  of  coordinates,  any  rotation 
of  the  input  images  cannot  be  accounted  for.  Not  only  would  the  degre- 
dation  of  the  maximum  peak  indicate  a  mis-match,  but  its  location  would 
also  be,  in  general,  a  biased  estimate  of  the  purely  lateral  displace¬ 
ment. 

Another  point  of  inconvenience  stems  from  the  large  amount  of 
computer  memory  and  the  inherently  long  CPU  time  required  to  evaluate 
the  inner  product.  The  calculation  of  a  single  inner  product  surface 
of  the  size  32  x  32  grid  points  from  images  of  the  size  128  x  128  takes 
on  the  PDP  11/60  typically  15-20  minutes  (pre-filtering  included). 
However  it  is  believed  that  a  much  faster  performance  could  be  achieved, 
if  special  function  processors  were  incorporated  in  the  computer 
hardware.  If  these  constraints  cannot  be  met,  new  displacement  esti¬ 
mators  have  to  be  developed.  These  estimators  then  should: 
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--identify  cloud  size  changes  quantitatively 
--identify  lateral  translation 
--identify  rotation 

--perform  fast  in  a  mini-computer  environment. 

The  estimator  introduced  in  this  chapter  is  based  on  fitting 
the  luminance  values  of  the  cloud  scene  (a  single  cloud  or  a  cloud 
field)  to  a  bivariate  normal  (Gaussian)  distribution.  The  proposed 
procedure  is  very  fast  and  yet  delivers  accurate  information  on  lateral 
displacement,  rotation,  and  size  changes. 


4.2  Theory 


The  luminance  or  brightness  values  (x(n,m)}  of  a  cloud  image  can 
be  fit  to  a  bivariate  normal  density  surface  placed  over  the  center  of 
luminance  gravity,  defined  as 


f(h,m)  =  B  exp< 


2(l-r2) 


(n-un)  ^  (fi-un)(m-um)  (m-uj 
o?  "  °n°m  o  2 


2  n 


(57) 


with 

1 

B  =  r - 7  ,  (58) 

27iao  i-r 
n  m 


Pn  and  ym,  the  center  of  luminance  gravity  (22),  and  the  second-order 
central  moments 


2  ,  N-l  M-l  ? 

=  ~  l  l  (n  -  u  )  x(n,m) 
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N-l  M-l 

E  E  (n  -  u  )(m  -  u  )  x(n,m), 
n=0  m=0  n  m 


(59) 
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with  Sqq  the  total  luminance  mass  of  the  image.  Note  that  f(n,m)  is 
defined  on  the  continuous  n,m  plane,  the  image  (x(n,m)}  however  is 
given  on  its  integral  grid  points.  The  term  of  (57)  enclosed  by  the 
square  brackets  indicates  that  the  locus  of  points  in  the  continuous 
n,m  plane  such  that  the  density  is  constant,  i.e.,  fvn,m)  =  c^,  is  an 
ellipse  of  the  normal  form 


(n  -  w„>‘ 


2r(n 


u  )  (m  -  u  ) 
n'v 


(m 


Mnr 


on, 

n  m 


a  2 
m 


'2* 


(60) 


with  its  center  (un»Vm)  as  1n  Fig.  6.  The  density  f(h,m)  is  maximum  at 
this  center,  and  it  equals 


f(Mn»Mm)  "  B,  (61) 

with  B  as  defined  in  (58). 

If  c2  =  1,  then  for  f(n,m)  =  f(un,m), 

(m  -  y  )2  =  a  2, 
m  m 

and  for  f(n,m)  =  f(n,ym). 


That  is,  the  radii  of  the  ellipse  parallel  to  the  coordinate  axes  are 
equal  to  on  and  om  (see  Fig.  6).  If  r  =  0,  then  the  principle  axes  of 
the  ellipse  are  parallel  to  the  coordinate  axes.  A  cloud  or  cloud  field 
is  fit  to  a  bivariate  normal  density  by  calculating  the  moments  and 
central  moments  of  the  image  brightness  and  equating  them  with  the  mo¬ 
ments  that  define  the  normal  distribution  of  luminance  values.  Changes 
of  the  parameters  of  the  ellipse  indicate  rotation,  lateral  displace¬ 
ment  and  size  changes  as  long  as  those  are  reflected  by  changes  of  the 
major  and  minor  axis  (denoted  by  a  and  b  in  Fig.  6). 


/  \ 

✓  \ 


'2 
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a  =  major  axis 
b  =  minor  axis 
a  =  angle  of  inclination 
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Figure  6.  Ellipse  of  constant  normal  density  for  c2  =  1. 


A  combination  of  this  process  and  the  inner  product  method  can 
be  developed,  where  the  rotation  is  estimated  first  by  identifying  the 
rotation  of  the  ellipse.  The  input  images  then  are  rotationally  re¬ 
aligned,  and  subsequently  the  inner  product  method  is  employed  for  a 
more  accurate  estimation  of  the  lateral  displacement. 

Once  the  moments  of  the  cloud  scene  are  established,  the  length 
of  the  principle  axes  and  the  angle  of  inclination  e  of  the  ellipse 
can  be  determined  as  follows: 

We  introduce  a  new  Cartesian  system  of  coordinates  placed  over 


the  center  of  the  ellipse. 
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with  parallel  to  the  major  axis  a,  and  parallel  to  the  rr.  1  nor  axis 
b  (see  Fig.  6).  IJithin  this  new  system,  the  coefficient  r  (see  (59)) 
should  be  equal  to  zero,  i.e.. 


_  I-  1-  4^2  x(t  ,c2)  =  0.  (63) 

r  =  —  »*  r  ~  —<xi 
'1  '2 

With  (62),  can  be  expressed  in  terms  of  n,tii  ana  (63)  becomes 

2  ?  2  2 
sin-  cost  u  +■  r(cos“  -  sin  ••):  >:  -  sin-  cost-  =  0. 

n  '  Min:  m 

Hence, 

,  2  m  o„ 

-  4  tan  — 2 - ?  (64) 

p  -  o 

n  iii 


The  length  of  major  and  minor  axis  can  be  obtained  by  solving  the 
following  equations. 

At  first,  from  the  properties  of  the  radii  of  an  ellipse  it 

follows  that 

2,2  2  2 

a  +  o  -  on  +  ,;m  .  (65) 

Furthermore,  the  pre-factor  B  of  (57)  must  be  the  same  in  the 
system,  i.e., 
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or 


a  b  = 


o  c 

n  m 


1  -  r‘ 


(66) 


Now  (65)  and  (66)  are  two  equations  that  can  be  solved  for  a  and  b. 
Note  that 


a  +  b  •  0:  for  any  case 
a  -  b  >  0:  for  r  >  0 


(67) 


a  -  b  «•  0 


for  r  ■-  0. 
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4 . 3 _  Examples 

One  of  the  operational  differences  between  the  inner-product 
method  and  the  method  of  tracing  ellipses  of  constant  gray-value  dis¬ 
tribution  is  that  the  latter  requires  the  ‘man-in-the-loop1  to  determine 
the  same  cloud  or  cloud  field  on  both  input  images.  Whereas  the  maxi¬ 
mum  of  the  inner  product  can  be  used  as  a  criterion  for  the  matching 
condition,  the  bivariate  normal  fit  requires  this  information  a-priori. 
The  method  is  applied  best  in  an  interactive  type  of  operation.  The 
cloud  analyst  selects  a  cloud  of  interest.  From  the  evaluation  of  the 
histogram  he  determines  a  threshold  and  positions  a  frame  around  the 
target  cloud.  The  computer  program  then  is  envoked  in  order  to  obtain 
the  parameters  of  the  ellipse  for  the  gray  values  within  the  selected 
frame.  The  results  are  obtained  numerically  as  well  as  utilized  for 
an  overlay  on  top  of  the  cloud  image  (see  Fig.  7). 


Figure  7.  Ellipse  of  constant  gray-value  density  as  overlay  on  cloud 
image.  Elliptic  parameters  are  determined  for  gray  values  within  the 
interactively  selected  frame. 
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A  time  series  or  these  parameters  then  is  established  by  analyzing  a 
sequence  of  images.  An  evaluation  of  this  time  series  might  support 
motion  analysis  and  facilitate  the  prediction  of  the  elliptic  parame¬ 
ters. 

As  an  example  we  considered  the  cloud  pairs  of  Fig.  2.  The 
parameters  of  the  ellipses  are  listed  in  Table  4.  The  resulting  values 
for  the  lateral  displacement  (i.e.,  the  displacement  of  the  center  of 
luminance  gravity)  should  be  compared  with  Table  3.  As  one  might 
expect,  the  results  are  rather  close  for  the  case  of  the  cloud  with 
only  few  changes.  In  the  assessment  of  the  other  cases  one  should  keep 
in  mind  that  the  number  of  computations  of  the  bivariate-fit  method  is 
considerably  lower  than  for  the  inner-product  method.  With  the  PDP 
11/60  equipment  described  in  the  previous  chapter,  the  results  for  a 
128  x  128  image  are  obtained  in  typically  5  ...  10  seconds. 

4.4  Summary 

In  this  chapter  a  method  was  introduced  that  delivers  not  only 
estimates  of  the  lateral  displacement  vector,  but  also  of  rotation  and 
size  changes.  These  parameters  are  obtained  from  fitting  a  bivariate 
normal  surface  to  the  luminance  values  of  a  cloud  or  a  cloud  field. 

By  tracing  only  the  ellipses  of  constant  density,  motion  vectors  and 
distribution  changes  are  identified.  Compared  with  the  inner-product 
method,  this  technique  not  only  delivers  additional  information  on  the 
cloud,  but  can  also  be  applied  for  whole  cloud  fields,  in  particular  if 
the  geometric  lay-out  of  the  scene  is  somewhat  elliptic.  This  method 
is  also  much  faster  than  the  inner-product  technique  and  requires  con¬ 
siderably  less  memory  allocation.  However  it  requires  the  a-priori 
information  of  match  by  the  human  operator. 
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Table  4.  Elliptic  parameters  for  images  of  Fig.  2. 


a  b 

6  major  minor 

inclination  axis 

in  degrees  in  pixel  units 


Cloud  Pair  I: 

Fig.  2  (a) 

180.2 

314.6 

24.3 

10.8 

4.2 

Fig.  2  (b) 

178.0 

316.1 

24.0 

11.1 

4.1 

Difference 

-  1.8 

1.5 

-  0.3 

0.3 

-0.1 

Cloud  Pair  II : 

Fig.  2  (c) 

214.0 

185.8 

27.30 

12.2 

3.8 

Fig.  2  (d) 

221.6 

184.4 

24.34 

19.2 

4.6 

Difference 

7.6 

-  0.6 

-  2.96 

7.0 

0.8 

Cloud  Pair  III: 

Fig.  2  (e) 

325.5 

135.0 

-  4.64 

14.1 

10.0 

Fig.  2  (f) 

336.8 

137.1 

-  3.04 

20.1 

12.1 

Difference 

11.3 

2.1 

-  1.60 

6.1 

2.1 

center  of 
el  1  ipse 
(pixel  units) 
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FOURIER  DESCRIPTORS  AND  ELLIPTIC  APPROXIMATION 
OF  CLOUD  CONTOURS 

5.1  Introduction 

The  approaches  taken  in  the  previous  two  chapters  demonstrate 
the  fact  that  any  method  of  motion  analysis  is  based  on  a  particular 
model  relating  the  given  image  to  a  meaningful  description  of  the  pic¬ 
torial  information.  As  we  have  seen,  the  inner-product  method  is 
originally  based  on  the  distribution  of  the  signal  energy.  The  pre¬ 
filter,  however,  operates  as  an  edge-enhancer  and  thus  shifts  the 
emphasis  to  the  information  embedded  in  the  edges  and  contours. 

The  method  described  in  the  previous  chapter  is  based  on  the 
derivation  of  a  simple  geometric  planar  curve,  the  ellipse,  from  the 
distribution  of  luminance  values  and  a  subsequent  analysis  of  the 
parameters  of  this  ellipse. 

In  this  chapter  now  we  shall  depart  one  step  further  away  from 
the  gray  level  information  of  the  image  by  introducing  a  parametric 
description  of  closed  contours.  The  only  remaining  link  to  the  gray 
level  Information  consists  of  the  fact  that  the  contour  is  defined  as 
the  line  of  constant  gray  level.  In  the  case  of  visible  data,  the 
contour  could  be  defined  by  a  gray  level  slightly  higher  than  the  dark 
background,  if  the  whole  cloud  is  to  be  encompassed.  In  the  case  of 
IR  imagery,  the  pixel  luminance  represents  temperature  and  usually  can 
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be  associated  with  altitude.  Hence  a  contour  could  be  defined  at  a 
particular  temperature  of  interest. 

Many  times  it  is  possible  to  find  distinct  contours  even  in 
situations  where  the  clouds  are  not  clearly  separated.  This  feature 
could  be  employed  by  the  cloud  analyst  who  wants  to  focus  on  isolated 
phenomena  like  overshooting  cloud  tops.  After  having  analyzed  the 
histogram  of  that  particular  area,  he  selects  a  threshold  and  obtains 
a  contour  around  an  area  with  temperature  values  relatively  lower  than 
outside  the  boundary. 

Recent  research  [35  -  41]  provides  evidence  that  the  Fourier 
descriptor  (FD)  is  a  simple  and  precise  method  for  describing  and  com¬ 
paring  the  contours  of  closed,  planar  figures.  Common  to  all  cited 
publications  is  the  derivation  of  a  functional  description  of  the  given 
contour  and  its  subsequent  Fourier  series  (FS)  expansion.  Interpreta¬ 
tion  and  analysis  of  the  contour  then  is  performed  in  the  frequency 
domain. 

With  respect  to  the  functional  description  of  the  contour  in  the 
image  plane,  two  types  of  FO's  have  been  developed.  The  older  method 
[36,37]  is  based  on  a  normalized  angle  versus  arclength  function.  In 
spite  of  its  fast  matching  algorithm,  its  practicabil ity  is  limited  for 
several  reasons.  First,  due  to  the  normalization  process  in  the  image 
domain,  the  resulting  FD  does  not  contain  all  shape  information.  Second, 
the  inverse  transform  does  not  necessarily  result  in  a  closed  contour. 
Furthermore,  as  Richard  and  Hemami  have  pointed  out  [38],  this  type  of 
FD  is  very  sensitive  to  the  noisy  fluctuations  inherent  in  a  fuzzy 
boundary. 
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The  approach  taken  in  this  work  is  based  on  a  complex- valued 
contour  description  and  essentially  the  discrete  version  of  the  FD 
suggested  by  Granlund  [35]  and  later  refined  by  Persoon  and  Fu  [39]. 

The  contour  is  traced,  thereby  yielding  a  complex,  periodic  function 
of  time.  The  FD  expansion  of  this  function  retains  all  shape  information. 

Two  contours  are  analyzed  and  compared  by  means  of  their  FD  para¬ 
meters.  Minimization  of  the  Euclidean  norm  of  the  FD  parameters  of 
two  contours  results  in  estimates  for  lateral  displacement,  rotation  and 
scaling  that  match  the  contours  in  the  optimum  sense.  In  many  cases  of 
application,  the  contours  are  sufficiently  smooth  such  that  the  FS  con¬ 
verges  rapidly.  In  those  instances,  the  matching  procedure  can  be 
undertaken  with  few  low-order  coefficients  only. 

The  inverse  transformation  of  a  truncated  FD  leads  to  closed 
contours  that  are  a  smoothed  least-squares  approximation  of  the  original 
contour.  It  will  be  shown  that  the  inverse  transformation  of  the  two 
lowest-order  coefficients  results  in  an  ellipse  fitting  the  original 
contour  in  the  least-squares  sense. 

It  should  be  noted  that  the  FD  is  a  method  to  describe  and 
analyze  general,  closed  contours  and  therefore  can  be  applied  not  only 
for  gray-level  contours  of  satellite  images,  but  also  for  isobars  and 
other  lines  connecting  points  of  equal  value  of  a  physical  parameter. 

5.2  Theory  of  the  Complex  Fourier  Descriptor 


5.2.1  Definition 

We  are  given  a  continuous,  closed  contour  C  in  the  continuous 
x-y-plane.  By  tracing  this  contour  with  constant  velocity,  a  complex, 
periodic  time  function 
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z(t)  =  x ( t )  +  j  y(t)  (68) 

is  defined  for  all  contour  points.  The  velocity  of  traversing  the 
contour  is  chosen  such  that  the  time  required  to  complete  one  period 
is  T  =  2tt,  with  arbitrary  starting  point.  Hence  the  fundamental 
angular  frequency  of  z(t)  is  equal  to  1. 

The  Fourier  descriptor  (FD)  of  C  then  is  defined  as  the  complex 
Fourier  series  expansion  of  z(t),  i.e., 

oo 

z ( t )  =  Y.  c(n)  exp(jnt)  (69) 

n=-<" 

where  the  complex  coefficients  c(n)  are  given  by 

c(n)  =  ~ 

1  -T/ 

The  coefficients  c(n)  can  be  expressed  in  terms  of  their  real  and 
imaginary  part  as 

c(n)  =  a(n)  +  j  b(n)  (71) 

or  in  terms  of  amplitude  and  phase  as 

c(n)  =  p(n)  exp( j<j>(n) ) .  (72) 

The  FD  can  be  viewed  as  the  power  spectrum  of  the  contour 
fluctuations.  Convergence  of  (69)  is  assured  at  any  t  and  uniformly 
on  any  closed  interval,  whenever  the  Dirichlet  conditions  are  satis¬ 
fied,  i.e.,  for  each  finite,  closed  contour  with  a  finite  number  of 
discontinuities.  Smooth  contours  result  in  rapid  convergence  of  the 
Fourier  series  coefficients  c(n).  Since  z(t)  is  a  complex  function, 
in  general  c*(n)  f  c(-n),  with  the  asterisk  denoting  the  complex 
conjugate. 


r 


z(t)  exp(-jnt)  dt. 


(70) 
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contains  a  listing  of  these  operations  in  the  image  (time)  and  the  fre¬ 
quency  domain.  Note  that  the  amplitude  spectrum  of  the  FD  does  not 
change  with  any  of  the  geometric  operations  listed  in  Table  5,  with  the 
exception  of  the  zero-frequency  term  c(0)  indicating  a  lateral  dis¬ 
placement  of  the  contour. 


Table  5.  Geometric  operations  and  the  Fourier  descriptor. 


Geometric  Operation  Image  Domain  Frequency  Domain 

1)  original  contour  z(t)  {c(n)}7« 

2)  lateral  shift  by  w(t)  =  z(t)  +  B  d(0)  =  c(0)  +  B 

Xl’yi  W1thB  =  xj  +  jy1  d(n)  =  c(n)  for  n  t  0 

3)  rotation  around  w(t)  =  z(t)  exp  (j<j>)  d (0)  =  c{0) 

centroid 

d(n)  =  c(n)  exp  (j<|>) 
for  n  ^  0 

4)  scale  change  w(t)  =  s  z(t)  d(n)  =  s  c(n) 

5)  shift  of  starting  w(t)  =  z{t  -  a)  d(0)  =  c(0) 

point 

0  <  a  <  2tt  d(n)  *  c(n)  exp  (jna) 

for  n  f  0 
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A  shift  of  the  starting  point  along  the  contour  results  in  a  phase 
shift,  as  it  is  well  known  from  the  Fourier  theory  of  real,  periodic 
functions. 

The  area  enclosed  by  its  FD  can  be  given  in  terms  of  the 
coefficients  c(n)  as  [39] 

°°  ? 

S  -  -  Z  jc(n) y  mr.  (75) 

n=-°° 


5.2.3  Matching  Contours  by  Means  of  the  FD 

The  match  of  two  contours  C  and  D  with  their  respective  complex 
functions  z(t)  and  w(t)  can  be  expressed  by  means  of  their  FD  series 
{c ( n ) }  and  (d(n)}.  As  a  measure  of  closeness  between  these  two  FD 
series  we  choose  the  Euclidean  norm 


d(C,D) 


CO 

l 


|c(n)  -  d(n)|2 


'  n=-°° 

According  to  Parseval's  theorem. 


1. 

'2 


-  w(t)|2dt. 


(76) 


(77) 


i.e.,  minimization  of  the  Euclidean  norm  corresponds  to  the  minimiza¬ 
tion  of  the  mean-square  difference  of  the  two  contours  z(t)  and  w(t). 

Another  important  result  of  the  theory  of  Fourier  series  states, 
when  applied  to  contour  functions,  is  that  the  contour  obtained  from  a 
truncated  Fourier  series,  i.e.,  the  function 
M 

Z-M  =  E  exp  (78) 

’  n=-M 

approximates  the  original  contour  z(t)  in  the  least  squares  sense. 

With  this  result  in  mind,  minimization  of  the  Euclidean  norm  can  be 
performed  on  the  truncated  FD  series  with  the  finite  sum 
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d(C,D) 


M  7  '"2 

E  |c(n)  -  d(n)  \ 
n=-M 


(79) 


Since  our  interest  is  directed  towards  estimation  of  rotation  6. 
scale  factor  s  and  starting  point  shift  a  that  would  map  a  contour 
z(t)  into  a  second  contour  w(t)  while  obtaining  a  minimum  Euclidean 
norm,  the  sum 

M  2 

J(C,D)  =  E  |c(n)  -  s  d(n)  exp  (j(3  +  na))p  (80) 

n=-M 

W*0 

has  to  be  minimized  with  respect  to  the  parameters  s,  3,  and  a.  The 
coefficients  c(n)  and  d(n)  of  (80)  are  the  FD  of  the  contours  z(t)  and 
w(t),  respectively.  Note  that  the  zero-order  coefficients  c(0)  and 
d(0)  are  not  included  in  (80)  since  they  indicate  the  lateral  shift 
of  the  centroid. 

The  following  is  a  brief  summary  of  Persoon  and  Fu’s  scheme  to 
solve  this  minimization  problem  [39].  Suboptimal  procedures  were  sug¬ 
gested  in  the  same  paper  [39]  and  by  Wallace  and  Mitchell  [40],  but 
are  not  considered  at  this  point. 

Set  c*(n)  d(n)  =  o(n)  exp  (jty(n)).  Minimization  of  (80)  then  is 
equivalent  to  minimizing 

J(C,D)  =  Ec(n)  c*(n)  +  s2  z  d(n)  d*(n) 

-  2  s  E  o(n)  cos(^(n)  +  na  +  3). 

The  summation  E  is  taken  over  the  indices  from  -M  to  M,  n  f  0. 
The  partial  derivatives  of  this  expression  with  respect  to  s,  3,  and 
a  then  are  set  to  zero,  and  we  obtain 


I  a(n)  cos(iMn)  +  net  + 
E  d(n)  d*(n) 


(81) 
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1 


tanB  =  - 


loin 

sin 

+  nal 

lain 

cos 

+  na; 

(82) 


and 


Io(n)  n  sinMn)  +  ng 
Wn)  + 


tan  ®  La(n)  n  cos 


na 


(83) 


By  combining  (82)  with  (83)  we  obtain  an  equation  in  a 

f(a)  =  Ea(n)  sin(i//(n)  +  na)  l  no(n)  cos(i^(n)  +  na) 

-  5lo(n)  cos(ijj(n)  +  na)  t  na(n)  sin(i^(n)  +  na).  (84) 

In  order  to  determine  the  optimum  values  s,  8  and  a,  we  first  solve 
f (a)  =  0.  Note  that  f(a)  is  a  periodic  function  with  several  roots 

f(a^)  =  0.  For  each  possible  shift  of  starting  point  a^ ,  the  angle  of 

rotation  3  and  the  scale  factor  s  are  calculated  according  to  (82)  or 
(83),  and  (81),  and  the  Euclidean  norm  (80)  evaluated.  The  set  of  8. 
a,  and  s  is  optimum  when  the  absolute  minimum  of  (80)  is  achieved. 

Since  (84)  cannot  be  solved  analytically,  a  numerical  root 
finding  algorithm  has  to  be  employed.  Since  the  FD  of  smooth  contours 
converge  rapidly,  the  summation  in  (81)  to  (84)  is  performed  in  those 
cases  only  over  a  small  number  of  terms. 


5.2.4  Oiscrete  Fourier  Descriptor  and  Aliasing  Effects 

In  all  practical  cases  of  application,  a  contour  is  given  as  a 
sequence  of  discrete  points  in  the  x-y  plane.  It  is  therefore  of 
interest  to  study  the  relationship  between  the  continuous  and  discrete 
Fourier  descriptors. 

Assume  that  we  are  given  N  points  of  a  continuous  contour  z(t). 
These  points  should  be  spaced  such  that  the  time  interval  between  two 
points  is  2tt/N.  With  (69)  the  FD  for  this  discretized  contour  then  is 

CO 

z(k^q-)  =  E  c(n)  exp(jnk^-),  k=0,  1,  .  .  .  N-l.  (85) 
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Since  the  index  of  the  summation,  n,  can  be  written  modulo  N,  i.e., 

n  »  m  +  rN,  with  m=Q,l  .  .  .  N-l 

r  any  integer,  (86) 

the  exponential  term  of  the  right-hand  side  of  (85)  can  be  re-written 

as 

exp(jnk2ir/N)  =  exp(j(m+rN)k27r/N)  =  exp(jmk2iT/N) , 
and  Eq.  (85)  becomes 

N-l  oo 

z(k2tr/N)  =  E  E  c(m+rN)  exp(jk(m+rN)27T/N) 

m=0  r=-oo 

or 

N-l 

z(k)  =  E  exp  (jmk2Tr/N)  c(m),  k=0,l,  ...  N-l  (87) 

m=0 

with  the  aliased  coefficient 

OO 

c(m)  =  E  c(m+rN) .  (88) 

f=-ot> 

Also  note  that 

c(-m)  =  c(N-m) .  (89) 

Equation  (87)  is  the  discrete  Fourier  descriptor,  based  on  a  finite 
sum  of  the  aliased  FD  coefficients  c(m).  Note  that,  in  general,  the 
original  coefficients  c(n)  and  the  aliased  coefficients  S(n)  are  not 
identical.  However,  if  z(t)  is  trigonometric,  i.e., 

M 

z(t) =  E  c(n)  exp(jnt), 
n=-M 

and  if  N  >  2M  +  1,  then 

c(n)  *  2(n)  for  |n|  <  M 
c(n)  *  0  for  |n|  >  M. 

The  aliased,  discrete  FO  coefficients  can  be  obtained  from  the  discrete 

N-l 

contour  points  {z (k) }q  as  follows.  For  the  finite  number  of  contour 
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points,  the  integral  (70)  becomes 
1  N_1 

g(n)  =  E  z{m27r/N)  exp(-jnm27r/N). 
m=0 

Now,  (87)  is  inserted  into  (90)  yielding 
,  N-l  N-l 

g (n )  =  m  2  Z  c(k)  exp(jkm2ir/N)  exp(-jmn27r/N) 
"  m=0  k=0 


,  N-l  N-l 

=  rr  Z  c(k)  Z  exp(jm(k-n)2TT/N). 
k=0  m=0 


The  right  summation  of  (91) 


exp(.jN(k-n)2tr/N)  -  1 
exp(3(k-n)27r/N)  -  1 

and  therefore 

g ( k )  =  ^  c(k)  •  N 

The  discrete  FD  then  is  giv 


n  be  written  as 
0  for  k  t  n 
N  for  k  =  n, 

by  the  pair  of  equations 


1  H~l 

c(k)  =  rr  I  z(m)  exp(-jkm2Tr/N) 
N  m=0 

N-l 

z(k)  =  Z  c(m)  exp(jmk2n/N) 
m=0 


(90) 


(91) 


(92a) 

(92b) 


5.2.5  Truncation  of  the  Discrete  FD 

The  objective  of  this  section  is  to  show  that  the  inverse  trans¬ 
form  of  a  truncated  FD  is  a  least  squares  approximation  of  the  original 
discrete  contour,  i.e.,  that 
M 

z'(k)  =  Z  £(m)  exp(jmk2n/N) ,  for  2M  +  1  <  N  (93) 

m=-M 


is  a  least-squares  approximation  to  (92b),  i.e., 
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J 


N-l  ? 
I  |z(k)  -  z'(k)|2 
k=0 


(94) 


is  minimum  for  z"(k)  as  given  by  (93). 

Assume,  z'(k)  is  not  given  by  (93),  but  by  the  FD  sequence  (d(k)}, 
i  .e. , 


M 

z'{k)  =  E  d(k)  exp(jmk2TT/N) , 
m=-M 


or,  in  terms  of  real  and  imaginary  parts,  by 
M 

z"(k)  =  Z  exp(jmk2TT/N)  •  (a(m)  +  jb(m)). 
m=-M 

With  z(k)  =  x(k)  +  jy(k),  the  function 


N-i |r  M  o  , 

J(x,y,a,b)  =  E  |x(k)  -  E  (a(m)  cos(^nk)  -  b(m)  sin(^rmk))l 
k=0<L  m=-M  N  N  * 

M 

+  [y (k)  -  I  (a(m)  sin(^mk)  +  b(m)  cos(^mk))J2 


is  to  be  minimized  with  respect  to  all  a(k),  b(k).  Setting  the  partial 
derivatives  of  J(x,y,a,b)  with  respect  to  a ( k )  equal  to  zero,  one 
obtains  after  some  lengthy  algebraic  operations 

I  E  a(m)  [cos(^Tjmk)  cos(^nk)  +  sin(%nk)  sin(^5nk)l  = 
k=0  m=-M  L  n  "  n  N  J 


N-l  ~  o_ 

E  x(k)  cos(~nk)  +  y(k)  sin(^Jnk)  (95) 

k*0  N  " 

Note  that  the  right-hand  side  of  (95)  is  equal  to  N  S(n),  with  a(n) 

the  real  part  of  the  aliased  coefficient  2(n).  Hence 

M  N-l  - 

E  a(m)  e  cos(^5k(m-n))  =  3(n). 
m=-M  k=0 


I 

N 


(96) 


r 


60 

Since  the  second  summation  in  (96)  is  equal  to  N  for  m  =  n,  but  0 
otherwise,  it  follows  that 

a(n)  =  a(n)  for  -M  <  n  <_  M. 

A  similar  result  can  be  obtained  for  b(n),  namely,  b(n)  =  b(n) 
for  -M<n<M.  We  therefore  conclude  that  the  2M+1  nonzero  coefficients 
d (n ) ,  -M<n<M,  of  a  discrete  contour  with  N  points,  N>2M+1,  that 
approximates  the  original  contour  of  N  points  and  N  FD  coefficients 
in  the  least-squares  sense,  are  identical  to  the  2M+1  low-order  FD 
coefficients  c(n)  of  the  original  contour,  i.e., 
d(n)  =  c(n)  for  -M<n£M. 


5.2.6  Elliptic  Approximation  of  Closed  Contours 

N-l 

We  are  given  a  discrete,  closed  contour  {z(k)}Q  with  FD 
N- 1 

coefficients  {c(n)}Q  .  Assume  that  the  FD  coefficients  are  truncated 

such  that 


(d(n)} 


N-l  _ 


(97) 


c(n)  for  n  =  0 
0  otherwise. 

The  FD  series  expansion  on  (d(n)}  then  equals 
y(k)  =  c(0) , 

i.e.,  the  contour  has  degenerated  to  a  single  point,  namely,  the  cen¬ 
troid  of  the  original  contour.  As  „e  have  seen  before,  the  zero-order 
coefficient  indicates  lateral  displacement  of  the  contour,  if  it  is 
evaluated  for  two  consecutive  cloud  images. 

Now  assume  that  the  truncation  is  performed  on  two  coefficients, 


i.e. , 


{d( n ) }q~ 1  - 


c(n)  for  n=0,l. 
0  otherwise. 


(98) 


The  expansion  on  these  coefficients  then  equals 
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y(k)  =  c (0 )  +  c ( 1 }  exp(jk-jj), 

i.e.,  a  circle  centered  at  c(0)  and  with  radius  |c(l)|.  From  the 
results  of  the  foregoing  section  it  is  apparent  that  this  circle  is  a 
least-squares  approximation  to  the  original  contour  z(k).  Since  rota¬ 
tion  and  starting  point  are  irrelevant  for  the  description  of  a  circle, 
the  first  two  FD  coefficients  could  be  employed  to  estimate  lateral 
displacement  and  scale  changes  between  two  contours.  Note  that  y(k) 
results  in  a  closed  figure  with  equidistant  contour  points. 

The  truncation  according  to 

N  ,  c(n)  for  n=0,l,N-l;  i.e.,  for  n=-l,0,l 
{d(n )}"_1  =  (99) 

0  otherwi se 

leads  to  a  curve 

y(k)  =  c(0)  +  c(l)  exp(jk-)  +  c( -1 )  exp(-jk^).  (100) 

For  k=0,l, .  .  .  N-l ,  (100)  is  identical  to  the  equation  of  N  points 
lying  on  a  general  ellipse  in  the  complex  plane,  centered  at  c(0).  The 
parameters  of  this  ellipse,  i.e.,  major  and  minor  axis  and  angle  of 
inclination,  are  not  as  readily  available  as  for  the  circle.  The  equa¬ 
tion  of  the  general  ellipse  in  parameter  form,  centered  at  the  origin, 
is  given  by 

z(s)  =  ~  exp(je)  exp(j(s-Sg))  +  ~  exp(je)  exp(-j (s+sQ) )  (101) 

with  s  the  parameter  0<s<2tt,  a  the  length  of  the  major  axis,  b  the 
length  of  the  minor  axis,  0  the  angle  of  inclination,  and  Sq  the  shift 
of  the  starting  point,  0<Sq<27t . 

The  parameters  of  the  ellipse  are  therefore  determined  from  the 
F0  coefficients  by  the  following  equations: 
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|c(-l)l  •*§* 

arg  [c(l)]  =  0  +  kQ 

arg  [c(-l )]  =  6  -  kQ  (102) 

with  kQ  the  integral  shift  of  the  starting  point.  It  appears  that 
the  elliptic  approximation  by  truncation  of  the  FD  can  be  used  to 
estimate  lateral  displacement,  rotation,  scale  change  and  shift  of 
starting  point  of  two  consecutive  contours.  One  should  keep  in  mind, 
however,  that  these  estimates  are  based  on  a  very  smooth  approximation 
of  the  contour,  and  that,  in  general ,  better  estimates  are  obtained 
with  a  higher-order  truncation  and  the  minimization  of  the  Euclidean 
norm  as  discussed  in  Section  5.2.3.  Nevertheless  the  elliptic  approxi¬ 
mation  delivers  estimates  very  fast  and  might  be  satisfactory  in  situa¬ 
tions  where  the  original  contour  has  assumed  a  somewhat  elliptic  shape. 

Note  that  the  points  described  by  Eq.  (100)  are  points  of  a  con¬ 
tinuous  ellipse,  but  that  these  points  are  not  equidistantly  spaced. 

In  other  words,  the  FD  of  a  discretized  ellipse  does  not  result  in  only 
three  non-zero  coefficients  c(0),  c(-l)  and  c(l).  The  parameter  s  in 
(101)  is  not  the  pathlength.  Pathlength  p(s)  as  a  function  of  this 
parameter  s  is  given  by  the  elliptic  integral  of  the  second  kind 


The  total  circumference  p(s=2tt)  =  A  corresponds  to  t  =  2n  seconds 
tracing  time.  Tracing  time  t  and  completed  pathlength  are  therefore 
related  by 
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t's)  =  ^p(s), 

with  s  the  parameter  of  the  ellipse  (101).  The  inverse  transform  of 
a  truncated  FD  (99)  therefore  does  not  result  in  an  ellipse,  but  in  a 
set  of  points  being  a  subset  of  the  ellipse.  The  effect  that  these 
points  are  not  equidistantly  spaced  along  the  continuous  elliptic  curve, 
could  be  viewed  as  the  Gibb's  phenomenon  of  the  Fourier  descriptor. 

5.3  Results 

The  FD  was  introduced  so  far  for  closed  contours  in  the  complex 
x-y-plane  (see  (68)).  Since  our  images  are  defined  along  an  n,m  matrix 
scheme,  the  notation  requires  a  slight  modification.  From  now  on  we 
shall  define,  without  loss  of  generality,  a  contour  as  a  sequence  of 
points  in  the  continuous  n,m  plane  as 

z(k)  =  m(k)  +  jn(k) ,  k  =  0,1,.  .  .  N-l.  (104) 

This  n,m  plane  is  identical  to  the  one  defined  in  Chapter  4.  The 
image  pixels  are  located  on  its  integral  grid  points. 

The  extraction  of  contours  is  usually  perceived  as  a  form  of 
edge  detection  and  often  based  on  statistic  assumptions  of  image  and 
image  noise  [42,43].  Since  our  contours  are  defined  as  lines  of  con¬ 
stant  gray  level,  they  do  not  necessarily  coincide  with  physical  edges. 
Also  most  of  the  contour  followers  found  in  the  literature  do  not 
always  guarantee  the  generation  of  a  closed  contour.  Because  the 
images  on  hand  are  usually  received  with  a  very  high  signal-to-noise 
ration,  image  noise  is  not  a  problem. 

Therefore  a  very  simple  gradient  technique  was  chosen  that  had 
been  applied  successfully  to  radar  echos  [44].  With  this  method  all 
points  of  the  fi.tft  plane  are  identified  where  the  specified  threshold 
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level  is  intersected  by  the  linear  interpolation  of  two  neighboring 
pixels. 

The  whole  contour  analysis  program  runs  interactively.  The  cloud 
analyst  selects  features  of  interest  from  a  displayed  cloud  image  and 
determines  the  threshold  along  which  a  contour  is  to  be  drawn,  by 
either  interpreting  a  histogram,  or,  as  in  the  case  of  IR  imagery,  by 
his  interest  in  a  particular  temperature  level. 

The  cursor  then  is  moved  into  the  area  which  is  expected  to  be 
the  inside  of  the  cloud  contour,  and  the  contour  search  initiated. 

At  first  the  starting  point  of  the  contour  has  to  be  found.  The 
search  is  started  at  the  position  of  the  cursor  and  extended  towards 
the  positive  m  (horizontal)  direction,  until  a  pixel  value  is  found 
that  is  lower  than  the  given  threshold.  The  location  of  the  first 
contour  point  then  is  determined  by  a  linear  interpolation  between  the 
pixel  values  inside  and  outside  the  contour,  and  its  intersection  with 
the  threshold  plane.  The  search  for  further  points  then  is  continued 
as  follows.  Each  point  of  intersection  leads  the  contour  either  into 
the  inside  of  a  2  x  2  pixel  square,  or  out  of  it  (see  Fig.  8). 

x(n,m)  x(n,m+l) 


Figure  8.  2x2  pixel  square  and  contour  points  of  entry  and  exit. 

If  an  entry  point  is  found,  a  search  for  the  exit  point  is  per¬ 
formed  counter-clockwise  along  the  border  lines  of  that  square.  The 


exit  point  determines  a  new  entry  point  for  the  neighboring  square, 
and  so  the  search  is  continued. 

The  algorithm  works  very  fast  (ca.  10  seconds  for  a  contour  of 
2000  points),  and  the  generated  contour  is  displayed  as  an  overlay  on 
top  of  the  cloud  image.  If  the  operator  chooses  to  select  this  contour, 
all  points  are  stored  in  a  file.  Otherwise  the  operator  can  move  the 
cursor  to  a  different  location  and/or  select  a  different  threshold  in 
order  to  generate  a  different  contour. 

Storage  of  the  contour  points  is  performed  in  form  of  the  Re¬ 
floating  point  values  of  the  coordinates.  Another  possibility  is  to 
convert  the  real  coordinates  to  the  nearest-neighbor  integers  and  store 
the  contour  by  its  chain-code  [45].  In  any  case,  the  contour  has  to  be 
re-sampled  to  obtain  real-valued  contour  points  in  the  n,m  plane  such 
that  their  spacing  is  constant.  Note  that  if  two  contours  are  to  be 
compared  with  each  other,  they  must  have  the  same  number  of  contour 
points,  evenly  spaced  along  the  contour.  That  coincides  with  the  re¬ 
quirement  of  tracing  the  contours  with  a  constant  speed  such  that  each 
circumference  is  done  in  2tt  seconds. 

The  process  of  re-sampling  therefore  starts  with  the  determination 
of  the  total  circumference  and  the  spacing  period.  For  sake  of  sim¬ 
plicity,  the  circumference  was  evaluated  along  the  piecewise  linear 
contour  given  by  the  entry  and  exit  points.  The  location  of  the  equi- 
distantly  spaced  new  contour  points,  however,  was  calculated  by  a  two- 
dimensional  quadratic  interpolation  scheme.  Performance  tests  with 
simple  geometric  figures  showed  that  the  amount  of  noise  that  was 
generated,  had  hardly  any  impact  on  the  motion  analysis.  The  advantage 
of  this  re-sampling  scheme  consists  of  the  computation  speed  that  is 
achieved. 


Figure  9  shows  a  typical  result  of  a  piecewise  linear  cloud 
contour  and  the  re-sampled  contour  points.  The  Fourier  descriptor  then 
is  obtained  by  a  Fast  Fourier  Transform  (FFT)  technique,  ("nure  10 
depicts  the  amplitude  spectrum  of  the  FD  as  obtained  from  c.ie  contour 
of  Fig.  9. 


Figure  9.  Example  of  piece-wise  linear  cloud  contour  and  re-sampled 
version  (depicted  as  single  points). 

Figure  10  demonstrates  the  rapid  convergence  of  the  FD  coeffi¬ 
cients.  It  also  becomes  apparent  that  the  amplitude  spectrum  of  the  FD 
is  not  symmetric,  since  the  contour  function  is  complex- valued.  The 
process  of  truncation  essentially  sets  the  coefficients  around  N/2  to 
zero,  with  N  the  total  number  of  contour  points.  The  inverse  transform 
of  a  truncated  FD  series  generates  a  smoothed  version  of  the  original 
contour.  This  effect  is  shown  in  Fig.  11  for  the  truncation  of  the  FD 
of  contour  Fig.  9.  In  the  first  case,  only  the  coefficients  c(0), 
c ( 1 )  and  c ( -1 )  ®  c(63)  were  not  set  to  zero.  The  resulting  curve  is  an 


Frequency 


Figure  10.  Absolute  value  of  the  FD  coefficients  obtained  from  the 
contour  of  Fig.  9.  N  =  64. 


ellipse.  In  the  second  case,  also  the  coefficients  c (2 )  and  c(-2)  = 
c(62)  were  utilized.  This  should  demonstrate  that  few  complex  numbers 
can  well  describe  general,  smooth  contours. 


Figure  11.  Inverse  transform  of  truncated  FD  series;  case  1:  3  non¬ 
zero  coefficients,  case  2:  5  non-zero  coefficients. 
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Two  contours  are  analyzed  for  lateral  shift,  rotation  and  scale 
changes  by  applying  the  optimization  technique  described  in  Section 
5.2.3.  The  FD  series  of  both  contours  are  calculated  and  the  order  of 
the  coefficients  selected  that  are  to  be  subjected  to  the  subsequent 
minimization  of  the  Euclidean  norm.  The  roots  of  Eq.  (84)  are  obtained 
numerically  by  means  of  Mueller's  iteration  scheme  of  successive  bi¬ 
section  and  inverse  parabolic  interpolation.  This  scheme  was  made 
self-adaptive  such  that  the  restriction  on  the  error  tolerance  for  the 
output  were  stepwise  reduced,  after  a  certain  number  of  iterations  had 
not  succeeded  in  finding  any  root. 

The  Euclidean  norm  was  evaluated  for  each  root.  The  set  of 
parameters  for  rotation,  scale  change  and  starting  point  shift,  that 
resulted  in  an  absolute  minimum  of  the  norm,  was  returned  as  the  final 
set  of  estimates. 

Table  6  lists  some  results  of  a  contour  analysis  of  the  cloud 
pair  II  of  Fig.  2.  Since  the  lateral  displacement  is  given  by  the 
zero-order  coefficient,  it  is  the  same  for  all  three  cases  of  trunca¬ 
tion.  Truncation  was  performed  as  follows:  for  case  I,  c(-l),  c(0) 
nonzero.  This  is  the  elliptical  approximation.  For  case  II,  c(-3) 

.  .  .  c(3)  non-zero.  For  case  III,  c(-IO)  .  .  .  c(10)  nonzero.  Since 
cloud  pair  II  of  Fig.  2  consists  of  well-separated  clouds,  the  threshold 
level  for  the  contour  was  set  slightly  higher  than  the  background. 

Table  6.  Numerical  results  of  the  FD  motion  analysis  of  cloud  pair  II, 
Fig.  2. 


Geometric  Operation 

Type  of  Truncation 

I  II 

III 

Lateral  Displacement,  An,  Am 

-2.2  1.6 

-2.2  1.6 

KO 

iH 

C\l 

CM 

1 

Rotation  in  Degree 

-2.6 

-1.2 

-0.9 

Scale  Change 

1.1 

0.9 

0.9 
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Table  6  illustrates  that  the  results  of  truncation  type  II  and 
III  are  not  dramatically  different  for  sufficiently  smooth  contours. 

The  whole  analysis  was  performed  with  the  PDP  11/60  equipment  described 
earlier  in  this  work.  Computation  time  never  exceeded  15  seconds  in¬ 
cluding  I/O-operations,  for  one  pair  of  contours. 

5.4  Summary 

The  method  of  motion  analysis  presented  in  this  chapter  is  based 
on  the  information  contained  in  closed  contours.  These  contours  can 
be  defined  either  as  encompassing  a  whole  cloud  or  by  following  a 
selected  level  of  pixel  value.  The  contour  is  traced  with  constant 
velocity.  The  coordinate  values  at  equidistant  points  along  this  con¬ 
tour  are  combined  to  a  complex  number.  The  contour  therefore  has  be¬ 
come  a  periodic,  discrete  time  function  in  the  complex  plane. 

The  Fourier  descriptor  (FD)  is  defined  as  the  discrete  Fourier 
series  representation  of  the  contour.  Lateral  displacement,  rotation 
and  scale  change  are  simple  operations  in  the  time  and  frequency  domain. 

Smooth  contours  result  in  FD's  that  converge  to  zero  rapidly. 

That  means  that  few  low-order  FD  coefficients  contain  most  of  the  shape 
information  of  the  contour. 

Two  contours  can  be  compared  by  minimizing  the  Euclidean  norm 
of  two  truncated  FD  series.  This  results  in  estimates  for  translation, 
rotation  and  scaling  that  would  map  one  contour  into  the  other.  The 
inverse  of  a  truncated  FD  series  approximates  the  original  contour  in 
the  least-squares  sense.  For  the  case  of  preserving  only  the  two 
lowest  non-zero  order  coefficients,  the  resulting  contour  consists  of 
points  on  an  ellipse.  These  points,  however,  are  not  any  more  evenly 
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spaced.  An  interactive  computer  program  was  developed  that  obtains  the 
estimates  of  motion  parameters  for  two  consecutive  contours  in  less 
than  15  seconds. 


CHAPTER  6 


TIME  SERIES  ANALYSIS  AND  ADAPTIVE  FILTERING 

6.1  Introduction 

In  the  foregoing  chapters  attempts  were  made  to  define  and 
extract  parameters  of  motion  from  a  given  pair  of  cloud  images.  These 
parameters  are  related  to  lateral  translation,  rotation  and  size 
change. 

This  chapter  now  deals  with  the  time  series  analysis  of  the 
parameters  obtained  from  a  sequence  of  consecutive  cloud  images.  It 
is  hoped  that  the  parameter  estimates  obtained  from  past  observations 
can  be  fitted  to  a  mathematical  model,  upon  which  the  prediction  of 
future  cloud  motion  and  related  weather  conditions  could  be  based.  The 
search  for  such  a  model,  or,  more  precisely,  a  class  of  models,  is 
currently  considered  to  be  one  of  the  most  urgent  research  topics  in 
the  atmospheric  sciences. 

Within  the  framework  of  this  research  no  attempt  is  being  made 
to  discuss  the  vast  amount  of  models  that  have  been  developed  for  more 
or  less  general  atmospheric  conditions  and  processes. 

Most  of  these  models  are  based  on  the  observation  and  measure¬ 
ment  of  variables  and  parameters  that  describe  the  physics  of  the 
atmosphere:  air  pressure,  temperature,  radiation,  wind  vectors,  etc. 

As  an  alternative  approach  to  those  very  often  complex  models, 
this  chapter  is  concerned  only  with  the  results  of  motion  analysis 
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techniques  as  described  in  the  previous  parts  of  this  work.  It  is 
assumed  that  besides  those  motion  estimates  no  further  information  is 
available. 

The  problem  therefore  consists  of  the  following  tasks: 

1.  select  one  or  more  targets  in  the  cloud  scene, 

2.  select  the  technique  of  motion  parameter  estimation, 

3.  obtain  estimates  of  these  parameters  from  past  satellite 
observations , 

4.  identify  a  model  that  fits  these  estimates  according  to 
some  criterion  of  optimality, 

5.  predict  future  values  of  these  parameters  on  the  basis  of 
the  model  and  past  observations,  and 

6.  if  desirable,  update  your  model. 


Assume  a  point  target  is  moving  in  the  x-y  plane.  Its  position 
x  and  y  at  times  t  is  denoted  by  x(t)  and  y(t). 

Assume  furthermore  that  x(t)  and  y(t)  can  be  modeled  by  poly¬ 
nomials  in  t  of  order  m,  and  additive  noise,  i.e., 

X(t)  =  Xq  +  Xjt  +  x2t2  +  .  .  .  +  +  vx(t)  (105  a) 

y(t)  =  y0  +  yxt  +  y2t2  +  .  .  .  +  ymtm  +  vy(t)  (105  b) 

with  Xq ,  .  .  . ,  xm,y0,  .  .  .,  ym  the  coefficients  and  vx(t),  vy(t)  the 
noise  processes.  Equation  (105  a)  and  (105  b)  may  be  combined  to  the 
vector  equation 

z(t)  =  H(t)  c  +  v(t). 


v 


(106) 
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with 


z(t)  = 


x(t)  ' 
y(t) 


the  observation  vector. 


H(t)  = 


fl  t  t2  ....  tm 

0 

i 

1 

1 

o  | 

1 

1 

1 

1 

1 

1 

1  t  t7  .... 

the  matrix 


of  exponentials,  and 

£  =  (x0  •  •  •  xm  J  yQ  •  •  •  ym)T  the  vector  of 
unknown  parameters. 

An  estimate  of  the  unknown,  c.,  can  be  found  by  the  least  squares 
technique.  This  method  essentially  produces  an  optimal  fit  of  a 
chosen  functional  form  (in  this  case  a  polynomial)  to  a  given  set  of 
measurements.  With  m  measurements  x(t^),  y(t^),  i  =  1,  .  .  .,  m,  we 
construct  the  vector  equation 


2  •  H  c  +  v. 


or,  in  long  hand. 


z(t  x)' 

H  (tj)' 

v(t  1y 

z(t2) 

H(tz) 

v(t2) 

= 

c  + 

z(t  ) 

— '  m  j 

H(tm) 

•  (107) 


The  least-squares  estimate  then  is  obtained  by  solving  [46] 

c  =  (fiT  Hf1  HT  z.  (108) 

Note  that  for  the  evaluation  of  the  estimate,  all  measurements  have  to 
be  utilized  together  at  one  time.  The  matrix  (HT  H)'1,  also  called 
the  pseudo- Inverse  of  H,  is  different  for  each  number  of  measurements. 
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but  it  can  be  calculated  before  any  of  the  observation  are  obtained. 

As  with  most  other  curve-fitting  techniques,  the  least-squares  method 
does  not  allow  any  insight  into  the  underlying  physical  nature  of  the 
motion  process.  Such  an  attempt  is  made  in  the  following  section. 

6.3  Dynamical  Motion  Systems  and  Recursive  Filtering 

Less  than  by  a  static  description  of  the  output  trajectory, 
many  physical  motion  processes  are  described  rather  by  linear,  dynamic 
equations. 

A  linear,  stationary  dynamic  system  is  given  in  state  space 
notation  by  the  following  vector  equation: 

x(t)  =  A  x(t)  +  B  u(t)  (109) 

and 

z(t)  =  H  x(t)  +  v(t),  (110) 

where 

x  the  nxl  state  vector 

u  the  mxl  vector  of  forcing  functions,  or  noise 
z  the  pxl  vector  of  observations 
v  the  pxl  vector  of  measurement  noise 
A,  B,  H  coefficient  matrices. 

Example. 

Consider  the  mechanical  system  shown  in  Fig.  12.  Let  c  be  the 
spring  constant  and  k  the  damping  coefficient  between  the  mass  m  and 
the  floor.  Newton's  law,  m  a  =  F,  gives 
m  y(t)  =  -  k  y(t)  -  c  y  +  u, 
with  u  the  external  force  (input). 


Let  us  choose  y  =  Xj  and  y  =  x2  as  state  variables.  The  dynamical 
state-space  description  of  the  system  then  is  given  by 


X1 

y  «  [i  o]  1 

x2 


Figure  12.  Example  of  dynamic  motion  process. 

Discrete,  dynamical  systems  are  governed  by  difference  equations.  A 
stationary,  linear  discrete  dynamical  system  is  given  in  state  space 
notation  as 

x(k)  =  *  x(k-l)  +  w(k-l)  (111) 

z.(k)  =  H  x.(k)  +  v.(k)  (112) 

with  $  the  state  transition  matrix.  If  a  continuous  system  is  given, 
but  the  observations  are  taken  in  discrete  time  instances,  the  whole 
system  can  be  discretized  and  the  transition  matrix  4>  derived  from  the 
matrix  A  of  Eq.  (109). 


The  state  vectors  of  both  continuous  and  discrete  dynamical  systems  can 
be  estimated  recursively  by  means  of  the  Kalman  filter.  Since  our 
cloud  image  observations  are  obtained  at  discrete  times,  we  shall  sum¬ 
marize  the  Kalman-filter  equations  only  for  the  discrete  case.  For  the 
system  given  by  (111),  (112),  assume 
w(k)  =  N(0,  Q(k)  ) 
and  v(k)  =  N(0,  R(k)  ). 

With  the  initial  estimates 
x(0)  =  E[x(0)] 

and  the  error  covariance  ’  (113) 

P(0)  -  E[(x(0)  -  x(0))(x(0)  -  i(0))T], 
the  extrapolation  or  prediction  of  the  state  estimates  and  the  error 
covariance  matrix  is  performed  according  to 


x(k+l  |  k)  =  <J>x(k)  (114) 

and 

M(k+1)  =  4>  P(k)  J  +  Q(k).  (115) 

With  the  new  measurement  zjk+1)  arriving,  the  new  state  estimate  is 
given  by 

x(k+l)  =  x(k+l  |  k)  +  K(k)  [z(k+l)  -  H  x(k+l  j  k)],  (116) 

with  K(k)  the  Kalman  gain  matrix  given  as 

K(k)  =  M(k)  HT  [H  M(k)  HT  +  R(k)]"1.  (117) 

The  error  covariance  is  updated  according  to 

P(k+1)  =  [I  -  K(k)  H]  M(k+1) .  (118) 

The  recursion  continues  with  Eq.  (114). 


The  Kalman  filter  is  the  optimum  linear  state  estimator  against 
an  expected-squared  error  cost  function  even  without  the  assumption  of 
Gaussian  noise  distribution  [47]. 
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6.4  Kalman  Filter  Tracking  of  Moving  Objects 

Recursive  tracking  algorithms  based  on  the  Kalman  filter  have 
entered  the  literature  about  10  years  ago  and  gained  considerable 
attention  since  then.  With  respect  to  the  transition  matrix,  the  pro¬ 
posed  models  may  be  divided  into  two  groups. 

The  first  type  is  usually  applied  to  the  estimation  of  the 
track  of  non-maneuvering  re-entry  vehicles  [48,  49].  The  model  used 
in  this  approach  is  closely  based  on  the  physical  equation  of  motion 
and  includes  typically  drag,  gravity,  centrifugal  and  Coriolis  force 
vectors  [48].  Since  some  of  these  parameters  are  unknown,  the  state 
vector  equation  becomes  a  non-linear  estimation  problem  and  is  usually 
solved  via  an  extended  Kalman  filter  [48,  50]. 

It  appears  to  be  desirable,  to  introduce  a  similar  approach  to 
the  motion  of  clouds.  The  horizontal  equation  of  motion  of  an  air 
parcel  is  given  in  pressure  coordinates  as  [51] 

dv. 

__  =  _v$  -  f  k  x  v  -  a  »,  (119) 

with  V$  the  horizontal  pressure  gradient  force,  f  Jc  x  v  the  Coriolis 
force,  and  a  v  the  frictional  drag;  v^  is  the  velocity  vector,  <t>  the 
pressure  field,  f  the  Coriolis  parameter,  k  the  unit  vector  in  vertical 
direction,  and  a  the  drag  coefficient. 

Even  if  we  assumed  that  the  motion  of  clouds  coincides  with  the 
dynamic  motion  of  air  flow,  we  would  again  obtain  a  non-linear  estima¬ 
tion  problem  and  would  have  to  resort  to  an  extended  Kalman  filter.  A 
thorough  study  of  the  feasibility  and  ramification  of  such  an  approach 
is  recommended  for  further  research.  An  extended  Kalman  filter  is 
likely  to  work  satisfactorily  for  a  very  large  number  of  observations. 
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However,  our  case  studies  on  hand  never  exceeded  12  measurements  of 
cloud  displacement. 

The  second  type  of  Kalman  filter  for  tracking  purposes  is 

derived  from  a  Taylor  series  representation  of  the  transition  matrix. 

The  position  x(t=T)  of  the  object  can  be  given  by  the  Taylor  series 

expansion  in  terms  of  the  derivatives  of  x(.)  as 

2  3 

x(T)  =  x(0)  +  T  x(0)  +  j-  x(0)  +  ~  x(0)  +  .  .  .  HOT,  (120) 
with  HOT  denoting  the  higher  order  terms.  The  series  expansion  (120) 
exists  for  continuous  trajectories  x(t). 

Assume  for  example  that  the  motion  process  can  be  modeled  in 
such  a  way  that  the  second  and  all  higher  derivatives  are  equal  to 
zero,  i.e., 

x(T)  =  x(0)  +  T  x(0) 
y(T)  =  y(0)  +  T  y(0). 

With  T  the  time  lapse  between  two  observations,  the  displacement  of  the 
target  then  is  governed  by  the  state  space  equation 
x.(k)  =  4>  x(k-l) 
or 
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If  there  remains  some  doubt  about  the  truncation  of  Taylor  series,  the 
residual  can  be  treated  as  an  additional  noise  term  [52  -  56].  If  it 
is  assumed  that  the  target  has  the  ability  to  maneuver,  this  residual 
noise  term  might  become  highly  correlated.  In  these  cases  an  augmented 
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state  vector  is  formed  with  the  correlation  coefficient  as  additional 
element  in  the  transition  matrix  [52,  54]. 

Common  to  all  Kalman  filter  tracking  techniques  are  the  follow¬ 
ing  disadvantages: 

1.  If  the  noise  covariance  matrices  Q(k)  and  R(k)  are  not 
exactly  known,  the  filter  becomes  sub-optimal. 

2.  Even  if  the  initial  state  estimates  are  well  placed,  the 
recursion  needs,  depending  on  the  complexity  of  the  system 
and  the  signal -to-noise  ratio,  a  relatively  large  number  of 
observations  in  order  to  converge.  In  the  case  of  radar 
target  tracking,  several  hundred  measurements  are  required 
to  identify  the  object's  path.  The  more  maneuvering  cap¬ 
ability  the  object  possesses,  the  complexer  the  tracking 
system  has  to  be  designed,  and  measurements  are  required  at 
an  even  higher  rate. 

6.5  Results 

At  the  beginning  of  the  research  it  was  hoped  that  cloud  motion 
could  be  modeled  with  respect  to  lateral  translation,  rotation  and 
scaling.  Sequences  of  images  were  obtained  depicting  a  variety  of 
weather  situations.  Since  it  was  anticipated  that  only  an  extended 
series  of  observations  could  produce  motion  forecasts  within  some  con¬ 
fidence  range,  only  hourly  IR  images  were  chosen,  since  they  are  re¬ 
ceived  throughout  the  night. 

As  method  of  parameter  extraction  we  chose  the  contour  analysis 
technique  by  employing  the  Fourier  descriptor.  Target  scenes  were  se¬ 
lected  and  their  contour  changes  analyzed  by  a  computer  program 
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throughout  the  sequence  of  images.  The  derived  motion  parameters  were 
listed  on  the  line  printer  in  form  of  time  series  data. 

Several  factors  might  have  influenced  the  mostly  unsatisfactory 
outcome  of  these  studies.  At  first,  many  times  a  single  cloud  could  be 
picked  up  only  for  3  to  4  consecutive  images.  At  the  end  of  this 
relatively  short  life  span,  the  cloud  had  merged  or  dissipated.  Some¬ 
times  the  identity  of  a  cloud  with  respect  to  earlier  images  was 
unclear. 

When  a  contour  was  selected  along  a  certain  temperature  level, 
the  fluctuation  of  the  shape  was  often  considerable.  This  effect  may 
be  due  to  the  unfavorable  complexity  of  cloud  dynamics  in  the  cases 
under  study,  and/or  to  diurnal  changes  of  temperature. 

One  of  the  advantages  of  the  FD  method  is  the  possibility  to 
generate  a  future  contour  synthetically  from  the  present  contour  and 
predictions  of  rotation,  lateral  displacement  and  scaling.  However, 
in  almost  none  of  the  cases  a  confident  prediction  could  be  made. 

The  best  results  were  achieved  in  those  few  instances,  when  the 
motion  of  the  selected  features  was  less  chaotic,  i.e.,  when  the  fea¬ 
tures  propagated  in  a  relatively  smooth  path  and  participated  all  in  a 
large  system  of  air  flow.  In  those  cases  the  parameters  of  interest 
could  be  picked  up  for  one  cloud  feature,  and  as  soon  as  this  feature 
had  become  inconsistent,  a  new  feature  was  searched  for,  in  the  hope  to 
obtain  an  uninterrupted  motion  analysis  for  the  governing  atmospheric 
system. 

It  has  become  apparent  that  more  research  has  to  be  undertaken 
in  order  to  demonstrate  the  feasibility  of  the  motion  analysis  and 
prediction,  as  it  was  intended  for  this  work. 


One  of  the  future  objectives  will  be  to  obtain  more  frequent 
observations.  It  is  true  though  that  the  cases  considered  for  this 


Table  7.  Coverage  of  features  and  temperature  level  of  contour. 


Image  received  at 


Features  and  Temperature  in  Celsius 


18:03 

19:03 

20:03 

21:03 

22:03 

23:03 

00:04 

01:04 

02:04 

03:03 

04:03 

05:03 

06:03 

07:03 

08:03 

09:03 


♦ 

© 

'  -54° 

no  observation 


* 

© 

-33° 


The  numerical  results  of  the  motion  analysis  are  listed  in  Tables  8  - 
12.  For  each  listed  time  instance,  the  differences  of  the  contour 
parameters  with  respect  to  the  earlier  observation  are  entered.  These 
parameters  are:  lateral  displacement,  rotational  and  scale  changes  by 
evaluating  only  2  F0  coefficients,  and  the  same  by  utilizing  10  FD 
coefficients. 


84 


Table  8.  Motion  analysis  for  feature  No.  1. 


lateral  dis¬ 
placement  2  FD  Coefficients  10  FD  Coefficients 


Time 

An 

Am 

Rotation 

Scale  Change  Rotation 

Scale  Change 

19:03 

-2.23 

11.2 

-3.3 

0.97 

-4.31 

0.96 

20:03 

-4.5 

4.86 

0.7 

1.04 

1.85 

1.03 

21:03 

-1.23 

6.07 

-1.4 

0.99 

-2.35 

0.97 

22:03 

1.85 

4.01 

-1.72 

0.99 

-1.48 

0.98 

23:03 

-6.9 

9.6 

-0.25 

0.97 

0.47 

0.97 

00:04 

-6.1 

13.3 

-5.25 

1.16 

-3.80 

1.01 

01:04 

02:03 

3.6 

-6.4 

3.63 

0.82 

6.05 

0.83 

03:03 

-6.4 

5.4 

2.45 

1.15 

1.75 

1.07 

04:03 

0.5 

-14.7 

10.1 

1.04 

10.51 

1.00 

05:03 

-0.1 

5.8 

-2.04 

0.95 

2.37 

0.93 

06:03 

-2.7 

12.6 

1.15 

1.11 

-0.27 

1.02 

07:03 

08:03 


09:03 
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Table  9.  Motion  analysis  for  feature  No.  2. 


lateral  dis¬ 
placement  2  FD  Coefficients  10  FD  Coefficients 


Time 

An 

Ani 

Rotation 

Scale  Change 

Rotation 

Scale  Change 

19:03 

-18.9 

16.7 

-0.22 

0.96 

-0.01 

0.95 

20:03 

19.1 

-11.5 

-10.16 

1.7 

-10.24 

1.68 

21:03 

2.5 

0.0 

0.4 

1.75 

-0.9 

1.76 

22:03 

23:03 

00:04 

01:04  -  no  data  available 

02:03 

03:03 

04:03 

05:03 

06:03 

07:03 

08:03 


09:03 
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Table  10.  Motion  analysis  for  feature  No.  3. 


lateral  dis¬ 
placement  2  FD  Coefficients  10  FD  Coefficients 

Time  An  Am  Rotation  Scale  Change  Rotation  Scale  Change 


19:03 


20:03 

12.9 

-9.1 

-2.62 

1.27 

-3.11 

1.21 

-9.2 

8.8 

2.58 

1.00 

3.56 

0.94 

-6.8 

8.4 

5.55 

0.91 

3.03 

0.91 

6.7 

1.9 

-6.43 

1.08 

-7.26 

1.06 

-2.7 

5.4 

-0 . 4? 

1.02 

-0.32 

0.97 

no  data 

02:03 

-0.2 

4.8 

-4.04 

0.74 

-5.17 

0.73 

03:03 

13.7 

-14.4 

-0.82 

1.73 

-0.92 

1.73 

-3.0 

9.3 

-3.18 

1.25 

-3.60 

1.19 

05:03 

-0.1 

2.0 

0.10 

1.45 

2.53 

1.46 

06:03 

0.8 

1.5 

0.37 

0.94 

2.84 

0.92 

-12.7 

18.8 

32.81 

3.03 

32.64 

2.97 

08:03 

09:0 
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Table  11.  Motion  analysis  for  feature  No.  4. 

lateral  dis¬ 
placement  2  FD  Coefficients  10  FD  Coefficients 

Time  An  Affi  Rotation  Scale  Change  Rotation  Scale  Change 

19:03 

20:03 

21:03 

22:03 

23:03 

00:04 

01:04  - no  data  available - 

02:03 

03:03 


04:03 

-  7.0 

5.6 

-11.24 

0.58 

-12.3 

0.55 

05:03 

-18.9 

11.9 

18.36 

0.69 

16.6 

0.67 

06:03 

0.5 

2.83 

-9.2 

0.81 

-  8.81 

0.80 

07:03 

-21.3 

13.2 

1.27 

1.50 

2.45 

1.46 

08:03 

-22.5 

11.6 

3.99 

0.97 

4.18 

0.94 

09:03 

-12.3 

6.32 

-1.45 

1.50 

-2.16 

1.48 

> 
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Table  12.  Motion  analysis  for  feature  No.  5. 


lateral  dis¬ 
placement  2  FO  Coefficients  10  FD  Coefficients 

Time  An  Am  Rotation  Scale  Change  Rotation  Scale  Change 


6.6  Summary 

In  this  chapter  we  studied  the  problem  of  modelling  motion 
processes.  As  long  as  the  trajectory  is  considered  as  being  a  poly¬ 
nomial  curve,  the  least-squares  estimator  can  be  applied.  This  method 
requires  all  past  observations  at  one  time.  Most  physical  motion 
processes  are  of  dynamical  nature  and  governed  by  differential  or  dif¬ 
ference  equations.  With  the  introduction  of  state-space  variables, 
the  input-output  relation  can  be  described  in  form  of  a  vector  equation. 
The  model  of  the  motion  process  is  either  derived  from  physical  equa¬ 
tions  or  from  a  Taylor  series  expansion.  Estimation  can  be  performed 
recursively  with  the  Kalman  filter.  In  the  case  of  a  physical  motion 
system,  the  estimation  becomes  non-linear  and  an  extended  Kalman  filter 
has  to  be  employed.  If  the  system  is  derived  from  the  Taylor  series 
expansion,  certain  assumptions  are  made  on  the  maneuverability  of  the 
target.  The  higher  order  terms  of  the  series  expansion  might  be  re¬ 
placed  by  noise  terms,  but  maneuver  noise  is  usually  correlated,  and 
the  state- vector  has  to  be  augmented. 

Finally  attempts  were  made  to  apply  any  of  these  methods  to  a 
time  series  of  motion  parameters  extracted  from  a  sequence  of  cloud 
images.  These  attempts  were  not  successful  since  the  «  ''eric 
processes  were  in  all  cases  of  a  very  complicated  nature.  Smooth  tra¬ 
jectories  existed  only  occasionally  within  a  range  of  3  or  4  observa¬ 
tions.  The  experienced  problems  were  discussed  and  recommendations 
for  further  research  made. 
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CHAPTER  7 


CONCLUSIONS 

The  objective  of  this  dissertation  was  the  analysis  of  cloud 
motion  from  satellite  images.  The  problem  was  approached  from  a  sig¬ 
nal  analysis  point  of  view.  Almost  no  consideration  was  given  to  the 
physical  parameters  and  properties  of  the  atmosphere. 

Motion  was  defined  as  a  change  of  those  image  parameters  that 
describe  a  cloud  or  a  cloud  pattern  most  distinctively.  The  lowest 
level  of  such  a  descriptive  signal  analysis  comprises  lateral  displace¬ 
ment,  rotation  and  size  changes  of  the  cloud. 

Current  cloud  motion  analysis  techniques  are  based  on  the  inner 
product  or  cross-correlation  method  and  deliver  estimates  only  for 
lateral  displacement.  Before  introducing  new  methods  in  this  research 
that  extract  also  rotational  and  size  changes,  an  improvement  of  the 
inner  product  method  was  developed. 

According  to  the  inner  product  method  the  lateral  displacement 
estimate  is  obtained  by  locating  the  maximum  of  the  inner  product  sur¬ 
face  of  two  consecutive  images.  This  maximum,  however,  is  very  often 
not  clearly  identifiable,  nor  is  this  method  selective  with  respect  to 
any  of  the  cloud  signal  features.  The  inner  product  should  be  evaluated 
only  for  those  features  that  are  relevant  for  a  signal  description.  The 
separation  of  the  clouds  from  the  background  can  be  considered  as  a 
fundamental  selection.  This  dissertation  shows  that  a  considerable 
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improvement  of  the  performance  of  the  inner  product  method  can  be 
achieved  by  applying  spatial  filters  to  the  input  images  prior  to  the 
inner  product  operation. 

In  order  to  generate  sharp  maximum  peaks  of  the  inner  product 
surface,  these  filters  are  developed  along  some  theoretical  considera¬ 
tions  as  inverse  filters.  Since  these  filters  possess  a  high-pass 
characteristic,  the  subsequently  evaluated  inner  product  is  based  on 
high-frequency  image  features  such  as  boundaries,  edges  and  texture. 
These  features  are  very  relevant  for  a  description  of  a  cloud  image. 

The  maximum  inner  product  then  indicates  the  energy  of  those  high- 
frequency  image  features  that  have  translated  laterally  in  a  homogen¬ 
eous  fashion.  If  the  cloud  changes  are  severe,  this  translation  is  not 
necessarily  homogeneous,  and  the  inner  product  might  result  in  a  mul¬ 
titude  of  relative  maxima. 

In  these  instances  the  filter  characteristic  can  be  changed 
gradually  from  high-pass  over  high-emphasis  to  all-pass,  such  that  the 
energy  of  lower  frequencies  has  more  and  more  influence  on  the  inner 
product  operation. 

The  pre-filtering  technique  was  tested  with  a  variety  of  input 
images.  Results  were  very  satisfactory.  The  only  constraint  to  this 
method  (besides  delivering  estimates  only  for  lateral  displacement)  is 
given  by  the  computer  power.  The  evaluation  of  the  inverse  filter 
coefficients,  the  filtering  process  and  the  subsequent  evaluation  of 
the  inner  product  surface  poses  a  serious  problem  to  the  users  in  a 
mini -computer  environment. 

In  the  following  two  chapters  new  techniques  were  Introduced  that 
also  provide  estimates  on  rotational  and  size  changes.  The  first  one 
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can  be  considered  as  almost  the  opposite  to  the  inner  product  method; 
it  is  a  very  fast  technique,  identifies  rotation  and  size  changes,  but 
is  based  on  a  very  general  description  of  the  cloud  signal,  namely,  on 
the  fit  of  a  normal  surface  to  the  luminance  values  of  the  image. 

Lines  of  constant  density  are  ellipses  around  the  center  of  luminance 
gravity.  Size  changes  can  be  identified  as  absolute  as  well  as  by  the 
excentricity  of  the  ellipse.  Rotation  is  indicated  by  a  changing  angle 
of  inclination.  This  method  is  not  necessarily  very  precise  but  can 
be  used  for  a  fast  over-all  evaluation  of  the  cloud  motion.  It  also 
can  be  applied  to  a  cloud  field. 

The  restrictions  imposed  on  the  next  technique  are  more  strin¬ 
gent.  The  meaningful  cloud  feature  is  now  defined  as  a  closed  contour, 
be  it  the  actual  boundary  of  the  cloud  or  a  contour  along  any  level  of 
constant  luminance  or  temperature  (for  IR  images). 

As  long  as  the  signal -to-noise  ratio  of  the  images  is  suffi¬ 
ciently  high,  a  simple  gradient  technique  proved  to  be  a  fast  and 
effective  contour  extractor.  The  analysis  of  the  contour  and  its 
changes  is  performed  by  means  of  the  Fourier  descriptor  (FD).  The  FD 
is  defined  as  the  discrete  Fourier  series  representation  of  the  closed 
contour.  Lateral  displacement,  rotation  and  scale  change  are  simple 
operations  in  the  space  and  frequency  domain.  Smooth  contours  result 
in  FD's  that  converge  rapidly.  That  means  that  few  low-order  coeffi¬ 
cients  of  the  FD  contain  most  of  the  shape  information  of  the  contour. 

Two  contours  are  compared  by  minimizing  the  Euclidean  norm  of 
two  truncated  FD  series.  This  results  in  estimates  for  translation, 
rotation  and  scaling  that  would  map  one  contour  into  the  other. 

The  inverse  transform  of  a  truncated  FD  series  approximates  the 
original  contour  in  the  least-squares  sense.  The  inverse  of  the  two 
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lowest-order  FD  coefficients  results  in  an  ellipse.  As  with  the  bi¬ 
variate  normal  fit,  these  ellipses  can  indicate  lateral  translation, 
rotation  and  size  changes.  More  precise  estimates  are  obtained  from 
FD  series  that  are  truncated  at  a  higher  order.  Assume  that  we  are 
given  predictions  of  lateral  translation,  rotation  and  size  change. 

With  the  FD  method  then  a  predicted  contour  can  be  generated  synthet¬ 
ically  from  the  FD  series  of  t.he  most  recent  observation.  Thus  the 
cloud  analyst  has  the  option  to  display  a  predicted  contour  and  compare 
later  with  the  actual,  updated  observation. 

An  interactive  computer  program  was  developed  that  obtains  a 
motion  analysis  from  a  pair  of  contours  in  less  thar>  15  seconds  with 
the  PDP  11/60  available  for  this  research. 

In  the  last  chapter,  models  for  motion  processes  are  studied. 

Like  most  other  physical  motion  processes,  cloud  motion  is  considered 
to  be  a  dynamical  process  that  is  described  best  in  state-space  nota¬ 
tion.  The  model  is  either  derived  from  a  Taylor  series  expansion  or 
from  physical  equations.  The  Taylor  series  approach  is  more  general 
and  leads  to  a  simpler  statement  of  the  estimation  problem.  More  re¬ 
search  is  recommended  to  utilize  physical  equations,  like  the  horizon¬ 
tal  equation  of  motion  that  is  used  to  calculate  wind  vectors  from  a 
balance  of  the  physical  forces.  Estimation  of  the  state-space  variables 
can  be  performed  recursively  with  a  Kalman  filter.  In  the  case  of  a 
model  based  on  physical  equations,  the  estimation  problem  might  become 
non-linear,  and  an  extended  Kalman  filter  has  to  be  derived.  In  spite 
of  the  mathematical  and  computational  complications,  such  an  estimator 
could  converge  faster,  since  it  makes  use  of  the  physical  properties  of 
the  atmosphere. 


Several  case  studies  of  cloud  image  sequences  were  analyzed  with 
the  Fourier  descriptor  method.  The  contour  parameters  of  various 
image  features  were  extracted  and  listed.  Attempts  to  find  smooth 
trajectories  succeeded  only  in  very  limited  cases.  Several  reasons 
might  be  responsible  for  this  failure.  First,  the  images  were  given  at 
a  rate  of  only  one  image  per  hour.  Many  features  therefore  could  be 
observed  only  over  a  span  of  3  to  4  images.  Secondly,  the  amount  of 
local  cloud  physics  and  disturbances  was  in  all  cases  larger  than  ex¬ 
pected.  Since  the  shape  of  the  temperature  contours  of  clouds  is  very 
sensitive  to  local  effects,  the  author  recommends  to  derive  contours 
not  for  a  particular  pixel  value,  but  for  a  range  of  values.  Local 
effects  might  thereby  be  averaged  out.  It  is  true  though  that  all 
cases  considered  for  this  work  depicted  severe  weather  situations  with 
admittedly  high  complexity.  Much  better  results  are  expected  with  a 
higher  image  rate,  and  more  experience  has  to  be  gathered  through 
further  case  studies. 
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APPENDIX 


"PHASE  UNWRAPPING  IN  TWO  DIMENSIONS" 


AI.  INTRODUCTION 


In  several  areas  of  one-  and  two-dimensional  signal  processing, 
the  phase  of  the  Fourier  spectrum  has  become  a  subject  of  considerable 
interest  [A2] .  The  importance  of  phase  in  image  processing  was  demon¬ 
strated  in  recent  papers  [A3,A4].  The  two-dimensional  phase  of  an 
image  conveys  usually  more  information  to  the  human  observer  than  the 
amplitude.  This  can  be  explained  by  the  relationship  between  phase  and 
the  location  of  spectral  energy  within  the  image.  The  human  observer 
seems  to  gain  more  information  from  the  location  of  spectral  energy  than 
from  its  amplitude.  Images  obtained  from  the  inverse  transform  of  the 
phase  term  only  show  superior  resemblance  to  the  original  than  those 
obtained  only  from  the  magnitude.  This  result  recently  led  to  the 
development  of  an  image  coding  scheme  emphasizing  the  phase  [A4]. 

A  numerical  evaluation  of  the  phase,  however,  is  commonly  ham¬ 
pered  by  the  ambiguity  of  the  arctangent  function.  Phase  unwrapping, 
i.e.,  the  evaluation  of  a  continuous  phase  curve  or  phase  surface,  is 
not  only  desirable  to  visualize  phase  behavior,  but  is  also  a  necessary 
prerequisite  for  the  computation  of  the  complex  cepstrum  if  it  was 
chosen  to  follow  the  method  of  computing  the  complex  logarithm  of  the 
spectrum  [A1,A5].  The  existence  of  two-dimensional  cepstra  was  shown 
in  [A6]. 

The  unwrapping  algorithm  proposed  in  this  correspondence  is  a 
two-dimensional  extension  of  Tribolet's  method  [Al].  It  incorporates 
a  suggestion  by  Bonzanigo  [A7],  namely,  the  substitution  of  the 
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direct-method  discrete  Fourier  transform  (DFT)  by  the  Goertzel  algo¬ 
rithm.  Typical  image  spectra  suffer  from  severe  leakage  and  aliasing, 
resulting  in  a  downgraded  performance  of  the  algorithm,  however.  As  a 
remedy,  we  shall  suggest  some  simple  pre-processing  techniques. 

All.  THEORETICAL  BACKGROUND 

Let  {x[n,m]}  be  a  finite  image  sequence  and 
N-l  M-l 

X(u,v)  =  F{x[n,m]}  =  E  E  x[n,m]  exp(-j27r  (Al) 

n=0  m=0 

its  discrete  Fourier  transform  (DFT),  briefly  referred  to  as  the 
spectrum  of  x[.].  Without  loss  of  generality,  the  image  boundaries 
are  set  N  =  M.  Note  that  X(u,v)  is  a  continuous,  generally  complex 
function  and  can  be  written  as  either 

X(u,v)  =  XR(u,v)  +  j  Xj  (u,v)  (A2) 

or  as 

X(u,v)  =  XA(u,v)  exp  {j  arg  [X(u,v)]>  ,  (A3) 

with  XR(.)  the  real  part,  Xj(.)  the  imaginary  part,  XA(.)  the  amplitude 
and  arg[X(.)]  the  phase  of  the  spectrum.  The  evaluation  of  amplitude 
and  phase  as  the  magnitude 

|X(u,v)|  =  (Xr2(u,v)  +  X^tu.v))*4  (A4) 

and  the  principle  phase 

X.(u,v) 

ARG[X(u,v)]  =  arctan  y-^-— %  (A5) 

R  * 

leads  to  the  following  ambiguities:  (i)  the  amplitude  is  forced  to  be 
non-negative.  Hence,  whenever  XA(.)  crosses  the  complex  plane  from  the 
first  to  the  third  quadrant,  passing  through  zero,  or  vice  versa,  the 
phase  experiences  a  discontinuity  of  ±  n.  (ii)  ARG  [X(.)]»  being  the 
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principle  phase,  leads  to  "wrap-around"  effects  modulo  2tt.  Figure  A1 
visualizes  these  two  conditions.  A  fictitious  phasor,  i.e.,  the  locus 
of  the  complex  spectrum  in  the  complex  plane,  experiences  a  phase  dis¬ 
continuity  of  +  tt  (case  (i).  Fig.  Ala)  and  of  -2tt  (case  (ii).  Fig.  Alb). 
Note  that  our  definition  of  the  arctangent  function  coincides  with  the 
ATAN2(B,A)  function  in  ANSI-FORTRAN,  and  not  with  ATAN(X) .  The  dif¬ 
ferences  are  illustrated  in  Table  Al. 


Figure  Al.  Discontinuities  of  the  principle  phase:  a)  case  (i), 
b)  case  (ii). 


Table  Al.  ATAN2  and  ATAN  FORTRAN  functions. 


Phase  (degr.)  ATAN2  ATAN 


0 

0 

0 

45 

45 

45 

90 

90 

-90 

135 

135 

-45 

180 

-180 

0 

225 

-135 

45 

270 

-  90 

-90 

315 

-  45 

-45 

360 

0 

0 

405 

45 

45 

The  derivatives  of  a  two-dimensional  spectrum  can  be  obtained  as 


The  right-hand  side  of  (A12),  however,  is  nothing  but  the  negative 
value  of  the  n-coordinate  of  the  center  of  (brightness)  gravity  of  the 
image.  Equivalently,  the  center  of  gravity  in  m-direction  is  equal  to 
the  negative  phase  slope  in  v-direction. 


AIII .  TRIBOLET'S  PHASE  UNWRAPPING 


Trlbolet's  unwrapping  algorithm  [Al]  for  the  one-dimensional 
phase  is  based  on  the  fast  Fourier  transform  (FFT)  of  the  data,  and 
the  FFT  of  the  data  premultiplied  with  the  index,  giving  the  complex 
conjugate  derivative  of  the  spectrum.  With  the  initial  phase  value 
defined  as  zero,  the  principle  value  is  compared  with  an  unwrapped 
phase  estimate  obtained  by  a  trapezoidal  integration  of  the  phase  de¬ 
rivative.  Whenever  the  difference  is  close  to  an  integer  multiple  of 
7f,  say  k*7T,  then  this  value  of  k - 7r  is  added  to  the  principle  phase  giv¬ 
ing  the  unwrapped  phase.  However,  if  the  difference  is  larger  than  a 
given  threshold  and/or  if  the  difference  of  consecutive  estimated  phase 
values  is  much  different  from  the  average  phase  derivative,  no  decision 
is  made  and  a  subprogram  called. 

Within  this  subprogram,  the  values  of  the  principle  phase  and 
the  phase  derivative  are  calculated  in  between  the  FFT  frequencies, 
i.e.,  {u)k  =  (2tt/N)  •  k,  k=0,  1,  .  .  .  N-l},  N  being  the  record  length. 

With  subsequently  increasing  subdivisions  of  one  interval  of 
neighboring  FFT  frequencies,  more  and  more  values  are  calculated  until 
the  behavior  of  the  phase  is  consistent  with  its  estimate.  Since  for 
a  phase  which  frequently  changes  its  direction,  the  number  of  addi¬ 
tionally  needed  DFT  calculations  can  become  very  high,  Bonzanigo  [A7] 
suggested  the  use  of  the  Goertzel  algorithm  [A2]  (including  a  phase 
correction  term)  instead  of  the  direct  method,  thus  reducing  the  number 
of  real  multiplications  by  half. 


AIV.  IMPROVEMENT  OF  PERFORMANCE 


Two  effects  are  likely  to  limit  the  performance  of  Tribolet's 
algorithm:  aliasing  and  leakage  (also  called  Gibb's  phenomenon). 
Aliasing  results  from  under-sampling,  leakage  from  the  finite  extension 
of  the  data  record  and  any  discontinuities  within  the  data  record 
[A8,  p.  58], 

Leakage  can  be  explained  as  the  convolution  of  a  sinc-function, 
i.e.,  the  spectrum  of  the  record-limiting  window,  with  the  generalized 
spectrum  of  an  infinite  data  record.  Thus  leakage  becomes  very  severe 
for  an  abrupt  termination  of  the  signal. 

Leakage  leads  to  ripples  in  the  spectrum  in  between  the  Fourier 
frequencies.  To  illustrate,  consider  the  simple  function  in  Fig.  A2a, 
with  Fig.  A2b  its  principle  phase  and  Fig.  A2c  the  phase  derivative 
evaluated  at  the  FFT  frequencies  according  to  (A8).  At  first  glance, 
the  phase  seems  to  be  linear  (in  accordance  with  the  phase  one  would 
expect  from  an  infinite  record),  hence  the  phase  derivative  Fig.  A2c 
should  be  constant.  Figure  A2d,  however,  shows  the  phase  behavior  in 
between  FFT  frequencies.  The  ripples  appear  and  indicate  why  the 
phase  derivative  of  the  signal  Fig.  A2a,  evaluated  at  the  FFT  frequen¬ 
cies,  is  not  constant.  Note  that  the  leakage  effect  in  this  case 
arises  mainly  because  of  the  overall  non-zero  signal  component.  In 
order  to  follow  ripples  and  a  generally  steep  phase,  Tribolet's  algo¬ 
rithm  can  perform  satisfactorily  only  if  the  thresholds  are  set  small. 
This  again  leads  to  numerous  additional  OFT  calculations  at  intermediate 
frequencies,  thus  imposing  CPU  time  problems. 

Unfortunately  it  appears  that  leakage  and  aliasing  are  common 
problems  in  image  processing.  Note  that  image  signals  are  usually 
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non-negative  and  terminated  abruptly  at  the  boundaries.  To  some  degree 
almost  all  digital  images  are  under-sampled. 

Leakage  can  be  reduced  by  making  the  transition  at  the  image 
boundary  less  “abrupt."  This  could  be  done  by  (i)  subtracting  the 
minimum  pixel  value  from  the  total  image,  thus  "lowering"  the  overall 
image  surface,  or  (ii)  data  windowing.  Furthermore,  the  phase  slope 
average  can  be  reduced,  and  thus  the  frequency  of  wrap-around  occur¬ 
rences  reduced,  by  shifting  the  center  of  gravity  into  the  data  origin. 
This  way,  part  of  the  data  is  shifted  into  the  negative  (non-causal) 
half  plane.  As  an  example,  consider  the  signal  given  in  Fig.  A3(a). 
Figure  A3(b)  shows  the  principle  phase  and  Fig.  A3(c)  the  phase  deriv¬ 
ative  at  the  FFT  frequencies  without  pre-processing.  Figure  A3(d)  was 
obtained  by  subtracting  the  minimum  value  of  x(n)  from  the  whole  record 
and  shifting  the  origin  into  the  center  of  travity.  To  obtain  the  FFT 
spectrum,  the  signal  values  about  the  negative  axis  can  be  attached  to 
the  right  end  (Fig.  A3(e)),  as  if  an  identical  period  is  following. 

Care  has  to  be  taken  only  with  the  pre-multiplication  of  the  data  in 
order  to  evaluate  the  phase  derivative  (Eq.  (A6)  and  (A8)).  The  values 
attached  to  the  right  end  are  to  be  premultiplied  with  the  negative 
indices  according  to  Fig.  A3(d).  The  same  caution  is  necessary  for  the 
DFT  evaluation  at  the  intermediate  frequencies.  Whereas  the  FFT  is 
identical  to  the  Fourier  series  of  a  periodic  extension  of  the  data, 
the  DFT  at  intermediate  frequencies  is  not,  but  rather  the  spectrum  of 
a  zero  padded  version  of  Fig.  A3(d).  Hence  the  negative  indices  for  a 
part  of  the  signal  have  to  be  taken  into  account.  Figure  A3(f)  now 
depicts  the  principle  phase  of  Fig.  A3(d)  and  Fig.  A3(g)  the  phase  de¬ 
rivative  (at  the  FFT  frequencies).  The  reduction  in  the  number  of 


(d)  pre-processed  signal 


Figure  A3.  Phase  unwrapping  of  a  pre-processed  signal. 
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phase  wrap-arounds  becomes  apparent,  and  with  Tribolet's  algorithm, 
including  the  modifications  as  described  above,  the  unwrapped  phase 
Fig.  A3(h)  was  obtained.  Since  the  displacement  of  the  center  of 
gravity  into  the  origin  has  the  effect  of  extracting  a  linear  phase 
component,  the  resulting  phase  sequence  is  ready  for  the  evaluation  of 
the  complex  ceps tr urn. 

AV.  TWO-DIMENSIONAL  PHASE  UNWRAPPING 

Consider  an  image  sequence  {x[n,m]}  being  shifted  in  the  n,m- 
plane  such  that  the  origin  n=0,  m=0  has  minimum  distance  to  the  center 
of  gravity.  With  its  region  of  definition  -K<n<L,  -I<m<J,  the  DFT 
becomes 

L  J 

X(u,v)  =  E  E  x[n,m]  exp(-j2ir(^-  +  ^))  (A13) 

n=-K  m=-I 

Note  that  the  DFT  is  defined  for  all  real  frequency  numbers  u  and  v. 

For  the  set  of  integer  numbers  within  the  u,v  -  plane,  we  introduce  the 
notation  [u],  [v],  and  the  spectrum  at  these  points  is 

X[u,v]  =  z  Z  x[n,m]  exp( ^  ))  (A14) 

n=-K  m=-I 

We  shall  refer  to  Eq.  (A14)  as  the  FFT-spectrum.  Since 

«p(-J2it  Jfi)  =  exp(-0‘2n  MM)  (A15) 

the  FFT  spectrum  (not  the  DFT)  can  be  written  for  positive  summation 
indices  as 

N-l  M-l  r.il  mFul 

X[u,v]  =  E  l  x[n,m]  exp(-j2u (%-*■  +  "m  ))  (A16) 

n=0  m=0  " 
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with  x[n,m]  the  periodic  extension  of  x[n,m]  (see  Fig.  A4a  and  A4b). 

A  mixed  "DFT/FFT"  is  defined  for  those  frequency  numbers  with  only  one 
coordinate  being  integer,  i.e., 


m=J  both  for  x[n,m]  and  x[n,m] 

(b)  image  x[n,m] 


Figure  A4.  Original  image  and  periodic  extension. 
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L  J  r  i 

X(u,[v])  =  E  E  x[n,m]  exp(-j2TT(1^- +  -^)) 
n=-K  m=-I  " 

L  M“  1 

=  E  exp(-j2n^)  E  xjn.m]  exp(-j27T^)  (A17) 

n=-K  m=0 

The  sequence  xm[n,m]  is  a  periodic  extension  of  x[n,;itj  only  in  in¬ 
direction.  Similarly, 

J  N-i  r  , 

X([u],v)  =  E  exp(-j2Tr^)  E  xn[n,m]  exp(-j2?rrij-^-‘-)  (A18) 

m=-I  n=0 

with  xn[n,m]  the  periodic  extension  in  n-direction. 

It  is  not  necessary  to  generate  three  different  periodic  exten¬ 
sions.  Even  the  mixed  DFT/FFT  can  be  performed  on  the  data  x[n,m]  by 
re-arrangement  of  the  summation  terms.  This  is  explained  here  only  for 
the  spectrum  X(u,  [v]): 

Define  the  vectors 


and 


Then 


Y(n;[v])  =  E  x[n,m]  exp(-j2^p-) 
m=0 


(A19) 


M-l 


Ym(n;[v])  =  E  xm[n,m]  exp(-j2Tr^-) 
m=0 


(A20) 


Ym(n;[v])  = 


Y(n;[v])  for  0<n<l 

Y(L-n;[v]}  for  -K<n<-1.  (A21) 

A  A 

Hence  the  row-FFT  Y(n;[v])  is  performed  on  the  set  x[n,m],  and,  for  a 
column-FFT,  the  output  does  not  have  to  be  re-ordered.  For  a  column- 
DFT,  however,  the  summation  terms  have  to  be  re-arranged  according  to 
Eq.  (A21). 


The  two-dimensional  phase  unwrapping  technique  proposed  in  this  work  is 
designed  as  follows.  Only  one  vector  (i.e.,  row  or  column)  in  the 
frequency  domain  is  unwrapped  at  a  time.  To  obtain  the  initial  unwrapped 
phase  values,  the  phase  is  unwrapped  for  the  first  column  ([u]  =  0)  or 
row  ([v]  =  0).  With  these  values,  the  phase  is  unwrapped  row  by  row 
(or  column  by  column)  along  integer  frequency  numbers.  Thus,  whenever 
the  DFT  is  to  be  calculated  at  intermediate  frequencies,  a  mixed  mode 
DFT/FFT  can  be  applied. 

To  obtain  the  phase  derivatives,  the  spectrum  of  n-x[n,m]  and 
m-x[n,m]  is  needed.  Note  that  n,m  are  defined  as  -K<n<L  and  -I<m<J. 
Hence,  to  obtain  the  FFT,  the  image  is  to  be  extended  periodically, 
but  the  negative  pre-multiplication  has  to  be  preserved.  The  mixed 
DFT/FFT  of  n-x[n,m]  then  is  given  as: 

L  J 

Zn(u,[v])  =  I  l  n*x[n,m]  exp(-j2ir(^- +  ^)) 

n=-K  m=-I 

L  0  r  i 

=  Z  n-exp(-j2rr2~)  l  x[n,m]  exp(-j27rP-j|p) ) 

n=-K  m=- I 

L 

=  Z  n-exp(-j2^)  •  Ym(n;[v]).  (A22) 

n=-K 

Hence,  pre-multiplication  is  not  a  part  of  the  row  FFT,  but  of  the 
column  DFT.  A  similar  expression  can  be  obtained  for 

L  J 

Zm([u],v)  *  Z  Z  m-x[n,m]  exp(-j2ir(^  +  ^)).  (A23) 

'  n=-K  m=-I 

Since  the  derivative  is  needed  only  with  respect  to  the  direction  of 
unwrapping  along  an  integer  valued  frequency  coordinate,  the  evaluation 
of  Zn([u].v)  and  Z  (u,[v])  is  not  needed. 
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The  following  list  is  a  summary  of  all  steps  introduced  above  (see 
Fig.  A5).  It  was  chosen  to  unwrap  the  phase  row-wise,  with  the  initial 
values  obtained  from  unwrapping  the  first  column. 

A)  Shift  center  of  origin  of  n,m-plane  to  center  of  gravity 

and  obtain  x[n,m]  by  periodic  extension  according  to  Fig.  A4. 

B)  Unwrap  first  column  arg  X[u,l] 

1)  do  row-FFT  on  x[n,m]  giving  Y([n];[v])  vectors. 

Actually,  since  only  Y([n];[0])  is  of  interest,  the  FFT 
reduces  to  the  simple  summation 

Y([n];[0])  =  V  x[n,m].  (A24) 

m=0 

2)  To  obtain  the  principle  phase,  phase  derivative  etc., 
perform  the  (column)  FFT  on  Y([n];[0]),  giving  X([u],[0]), 

/V 

and  also  on  n  •  Y([n];[0]),  with  the  factor  n  according 
to  the  re-ordered  scheme,  giving  Zn([u],[0]). 

3)  If  the  (column)  OFT  is  needed,  i.e.,  X(u,[0])  and 
Zn(u,[0]),  a  re-ordered  summation  of  Y([n],[0]),  and 
n  •  Y([n];[0])  is  required. 

C)  Unwrap  all  phase  rows  arg  X[u,v] 

1)  column  FFT's  on  x[n,m] 

2)  do  row  by  row  [u] 

--look  up  starting  value 
-FFT  -  X[u,v] 

— pre-multiply  correctly,  and  FFT  +  Zm[u,v] 

--if  OFT  Is  needed:  re-order  m-summation,  and  OFT's 
give  X([u],v)  and  Zm([u],v). 
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Since  unwrapping  is  perform'd  vector-wise,  Tribolet's  algorithm  can 
be  utilized  with  the  following  minor  changes: 

--the  input  data  is  complex 

--pre-multiplication  and  DFT  summation  terms  need  occasional 
re-ordering  of  indices. 

AVI.  RESULTS 

By  shifting  the  image  origin  to  the  nearest  pixel  location  of 
the  center  of  luminance  gravity,  phase  plots  were  obtained  with  a 
derivative  of  almost  zero  at  the  spectral  origin. 

During  the  generation  of  various  phase  surfaces,  however,  it 
was  discovered  that  the  smoothness  of  the  surface  plot  could  be  dras¬ 
tically  improved,  and  thus  the  number  of  DFT  calculations  minimized, 
by  choosing  a  different  point  of  the  image  as  the  spatial  origin. 

This  point  is  to  be  selected  such  that  the  unwrapped  phase  at 
the  end  points  of  the  frequency  range  are  close  to  zero.  With  the 
signal  and  spectrum  related  to  each  other  by  Eq.  (Al),  a  shift  of  the 
image  by  n^,  mQ  is  equivalent  to  a  linear  phase  shift,  according  to 

nnu  mnv 

x(n-n0,m-m0)  =  X(u,v)  exp(-j2ir(-^- +  -jj-)).  (A25) 

The  required  phase  shift  then  is  determined  by  evaluating  the  unwrapped 
phase  arg  [X(Q,M)]  and  arg[X(N,Q)].  Since  these  values  are  calculated 
by  essentially  one-dimensional  unwrapping  procedures,  they  are  obtained 
fast. 

In  addition  to  this  improvement  it  was  desired  to  obtain  a  phase 
surface  with  its  spectral  origin  at  the  center.  This  was  achieved  by 
exploiting  the  conjugate  symmetry  of  the  spectrum. 
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Figure  A6  serves  as  an  example  for  such  an  unwrapped  phase  sur¬ 
face.  Note  that  the  unwrapped  phase  is  very  sensitive  to  estimation 
errors  in  spectral  regions  of  low  energy,  i.e.,  where  the  phasor  is 
close  to  its  origin  [A9]. 


Figure  A6.  Example  of  unwrapped  two-dimensional  phase. 
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